A brief compendium of time-dependent density-functional theory 
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Time-dependent density-functional theory (TDDFT) is a formally exact approach to the time- 
dependent electronic many-body problem which is widely used for calculating excitation energies. 
We present a survey of the fundamental framework, practical aspects, and applications of TDDFT. 
This paper is mainly intended for non-experts (students or researchers in other areas) who would 
like to learn about the present state of TDDFT without going too deeply into formal details. 



I. PREFACE 

This paper presents an introduction to and a survey 
of time-dependent density- functional theory (TDDFT). 
The purpose of the paper is to explain in a nutshell what 
TDDFT is and what it can do. We will discuss the basics 
of the formal framework of TDDFT as well as the cur- 
rent state of the art, skipping over details of the proofs, 
and highlight some of the most important applications. 
Readers who would like a more detailed treatment and 
more literature references are encouraged to consult re- 
cent books [H Q and review articles [Jll] ■ 

TDDFT is a theoretical approach to the dynamical 
quantum many-body problem; it can be used to describe 
quantum systems that are not stationary. As a conse- 
quence, TDDFT provides formally exact and practically 
convenient methods to calculate electronic excitation en- 
ergies. By contrast, density-functional theory (DFT) is a 
ground-state theory: in other words, it is used to find the 
ground state of a quantum system and calculate related 
quantities of interest, such as the ground-state energy. 
In many, if not most, situations of practical interest, we 
have to determine the ground state of the system before 
we can study its dynamics or calculate its excitations. 

The beginnings of ground-state DFT date back to the 
years 1964/65 when the famous papers by Hohenberg and 
Kohn ^] and Kohn and Sham [7,] were published. Since 
then, DFT has developed into a dominating method for 
electronic structure calculations in physics, chemistry, 
materials science, and many other areas (see Ref. i8] for 
a recent up-to-date account of DFT). Although TDDFT 
is of much more recent origin [Oj], it now has reached a 
similar status for calculating electronic excitations. 

TDDFT uses many familiar concepts from DFT, most 
prominently, the Kohn- Sham idea of replacing the real 
interacting many-body system by a noninteracting sys- 
tem that reproduces the same density. But there are also 
many concepts that are unique to the time-dependent 
case, such as memory and initial-state dependence. To 
gain a thorough understanding of TDDFT it is hence 
advisable to begin with a study of the basic concepts of 
DFT. We refer the reader to the very nice introductions 
to DFT by Capelle O and by Burke and Wagner IH. 
There exist a number of books on DFT, some of which are 
very accessible to newcomers in the field [l^, E^l ; others 
are more advanced [14| . 



II. GROUND-STATE DFT IN A NUTSHELL 

A. The many-body problem 

We consider a system of N interacting electrons that 
is described by the nonrelativistic Schrodinger equation: 

Ha-^j{vi,...,VN) ^ Ej-^j{vi,...,vn) , j =0,1,2,... . 

(1) 

For a given _D-dimensional system, the jth eigenstate 
'i>j{ri, . . . , r^r) is a function of DN spatial variables. In 
the following, we use the abbreviation Vfj. Of particular 
interest to us is the ground-state wave function ^fgg. 
The many-body Hamiltonian is given by 



H = T + V + W . 



(2) 



where the kinetic-energy and scalar potential operators 
are 



TV 



N 



2 ' 



(3) 



and the electron-electron interaction operator is 



N 



W=J2 ^(Ir.-r,|) 



(4) 



Notice that we use atomic (Hartree) units throughout, 
i.e., m = e = h = 1. The electron-electron interaction is 
usually taken to be the Coulomb interaction, i(;(|r— r'|) = 
l/|r — r'l, but other forms of two-particle interactions, or 
zero interactions, are also possible. 

The single-particle potential v{r) describes the total 
potential acting on the electrons. If one is interested in 
describing the properties of matter (atoms, molecules, 
or solids), v{r) is the sum of the Coulomb potentials of 
the atomic nuclei. However, to define the formal frame- 
work of DFT it is not necessary to specify where the 
potential comes from, as long as it has a mathematically 
well-behaved form. 

From the solutions of the Schrodinger equation we cal- 
culate the expectation value of an observable in the jth 
eigenstate: 



(5) 
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Here, O is a Hermitian operator corresponding to a quan- 
tum mechanical observable. 

Let us make two remarks on our formulation of the 
many-body problem. 

(a) We implicitly made the Born-Oppenheimer approx- 
imation (see Section fVlI Cp . In other words, if our system 
contains nuclear degrees of freedom (as is the case in all 
forms of real matter), we treat them classically. The 
many-body wave functions therefore depend only on the 
electronic coordinates (ri, . . . , r^v), and the nuclei only 
act as sources of scalar potentials. In Section I Villi we 
will briefly discuss what happens if this approximation 
is not made and the electronic and nuclear degrees of 
freedom are coupled. 

(b) We have not indicated any spin indices, which was 
done mainly for notational simplicity. In other words, 
^'j(ri, . . . , Tat) describes spinless electrons. Including 
spin, the many-body wave function can be written as 
^'(xi, . . . , xat), where x^ = {ri,ai) denotes the spatial 
and spin coordinate of the ith electron. 



B. The basic idea behind DFT 

Everything we wish to know about our system (en- 
ergy, geometry, excitation spectrum, etc.) can be ob- 
tained from the wave functions, see Eq. ([5|). The exact 
wave functions can be calculated if the system is small, 
with no more than one or two electrons, but this becomes 
very difficult if N is greater than that: the many-body 
Schrodinger equation becomes too hard to solve, and the 
usefulness of the wave function itself becomes more and 



more questionable for large N |15|. 

There exist many approaches to find approximate so- 
lutions of the many-body problem. So-called "wave- 
function based techniques" such as Hartree-Fock (HP) or 
configuration-interaction (CI) attempt to find variational 
solutions of the Schrodinger equation using expansions of 
the wave function in terms of Slater determinants. This 
approach has been very successful in theoretical chem- 
istry, but has its limitations for large systems. 

The essence of the density-functional approach is that 
it is in principle possible to obtain all desired information 
about an N -electron system without having to calculate 
its full wave function: instead, all one needs is the one- 
particle probability density of the ground state, 



no(r) =N / dVa... / dVjv|*gs(r, rj, . . . , r^) 



(6) 



This can be mathematically proven (see below), but 
before doing so it is helpful to give a simple illustra- 
tion. Consider a one-electron system which satisfies the 
Schrodinger equation 



' 2 



(7) 



The usual procedure is to solve this equation for a given 
potential v(r) and determine the ground-state probabil- 
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FIG. 1. The density no{x) — cos^(7ra::) -(- cos^(37r2;) (dashed 
line, scaled by a factor 60), for — | < a; < |, is the ground- 
state density of the potential v{x) (full line) which was con- 
structed using Eq. (|8|. 



ity density as no(r) — |(po(r)p. But now imagine the 
reverse situation: we are given a density function rio(r), 
normalized to 1, and we ask in what potential this is the 
ground-state density. Assu ming that the wave functions 
are real so (po (r) = •y/no(r), the Schrodinger equation ([7]) 
is easily inverted, and we obtain 



v{r) = 



_ V2no(r) |Vno(r) 



4^0 (r) 8no(r) 



(8) 



A one-dimensional example is given in Fig. [TJ 

What has been accomplished? From the ground-state 
density uqIj) we were able to reconstruct the potential 
v(r). But this means that we have reconstructed the 
Hamiltonian H of the system, and we can thus solve the 
Schrodinger equation ([7]) and get all the wave functions! 
This logical chain can be represented as follows: 



no(r) -> v{r) H ^ 



(9) 



The reconstruction of the potential from the density is 
easy for one-electron systems. For interacting systems 
with many electrons there is no explicit formula such as 
Eq. ([8]). Nevertheless, there exists a unique potential for 
each mathematically well-behaved density function such 
that it is the ground-state density in this potential. This 
was proved by Hohenberg and Kohn in 1964 Q. 

The Hohenberg- Kohn theorem states that it is impossi- 
ble for two different potentials, w(r) and ^'(r), to produce 
the same ground-state density [v' is considered to be dif- 
ferent from V if it is not just v shifted by a constant). 
In other words, the relationship between potentials and 
ground-state densities is one-to-one: 



no(r) o t;(r) 



(10) 



The proof of this theorem is relatively straightforward, 
making use of the Rayleigh-Ritz minimum principle. It 
can be found in any textbook on DFT, so we won't repeat 
it here. 

Formally, this logical dependence of the wave functions 
on the ground-state density constitutes a functional rela- 
tionship, which is written as 4'j[no]. Hence, the name 
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density-functional theory. Every quantum mechanical 
observable thus can be written as a density functional. 

In particular, the total energy functional of a system 
with potential vo(j) is 

E^,[n]^ {^[n]\f + Vo + W\^[n]) , (11) 

where n is some iV-electron density and ^[rt] is that 
ground-state wave function which reproduces this den- 
sity. The energy functional (llOp is minimized by the 
ground-state density uq which belongs to vq, and then 
becomes equal to the ground-state energy: 

E^^[n] > Eo for n{r) ^ no(r) , 

^,,[n]=So for n(r)=no(r). (12) 

C. The Kohn-Sham approach 

The fact that all observables are functionals of the den- 
sity opens up the way for an enormous computational 
simplification, since the density is a function of only D 
variables (and not of DN variables as the wave function). 
But how can we take advantage of this in practice? To 
obtain the density one still needs to solve the full many- 
body problem. This means that nothing has been gained, 
unless we find a way to bypass the full Schrodinger equa- 
tion and obtain the density in some other, easier way, at 
least approximately. Fortunately, a very elegant method 
exists to do this, known as the Kohn-Sham formalism [t']. 

We use the following trick: we define a noninteracting 
system in such a way that it reproduces the exact ground- 
state density of the interacting system. This means that 
we can calculate the exact density as the sum of squares 
of single-particle orbitals. 



N 

no(r)=^|^,(r)p, 
i=i 

where the orbitals satisfy the following equation; 
2 



- Vs [n]{r) 



(13) 



(14) 



This equation is known as Kohn-Sham equation; it is for- 
mally a single-particle Schrodinger equation, like Eq. ([7|) . 
However, the potential Vg is very special: it is defined to 
be that single-particle potential that produces orbitals 
which give the exact ground-state density of the inter- 
acting system via Eq. It is therefore a functional 
of the density, Vg [n] (r) . 

The trick is now to write the unknown effective poten- 
tial Vs [n] in a smart way. No doubt, the given external po- 
tential vo{r) will make a contribution to it. The remain- 
der, Vs[n\ — vo{r), then accounts for the electronic many- 
body effects. A large portion of the latter is made by 
the classical Coulomb potential associated with a given 
density distribution, also known as the Hartree potential. 



WH(r) 



-3 / "(!■') 



(15) 



And whatever is left is called the exchange- correlation 
(xc) potential, Vxc ["■](!■), so that 



Vs[n]{r) = vo{r) + WH(r) + WxcM(r) 



(16) 



It turns out that the solution of the Kohn-Sham equation 
(fT4|) — that is, the density (fT3|) — is precisely that density 
which minimizes the total energy functional ()lip . The 
connection is made by rewriting the energy as follows: 

Ey„[n] = T[n] + J d^r t;o(r)n(r) W[n] 

= Ts[n] + j d\vf,{v)n{v) + (T[n] - T,[7i] + W[n]) 

= r,[n] / dVvo(r)n(r) + EhM + E^,[n] . (17) 



Here, T[n] is the kinetic-energy functional of an inter- 
acting system, whereas Ts[n] is the kinetic-energy func- 
tional of a noninteracting system. Neither T[ n\ nor Tg \n\ 
are known as explicit density functionals, but it is very 
straightforward to write down Ts [n] as an explicit func- 
tional of the orbitals: 



1 ^ 



(18) 



where the orbitals fjijc) come from the Kohn-Sham equa- 
tion (|14p are are hence implicit density functionals. 
In the last line of Eq. PT)) we define the Hartree energy 



I n{r)n{r') 



(19) 



and the xc energy as 

[n] =T[n]- Ts [n] + W[n]- E^ [n] . (20) 

This shows that the xc potential is given by the following 
functional derivative: 



t'xc(r) = 



5E^ 



5n{v) 



(21) 



It is straightforward to show that the total energy p7)) 
can be expressed as follows: 



N 

Evo [n] = X! ^ + / '^^^ Wxc(r)n(r) + E^c[n] 



i=i 



D. Discussion and exact properties 



(22) 



Let us now summarize some of the most important 
properties of the Kohn-Sham approach. Our discussion 
is by no means complete, but the following properties will 
be relevant for the time-dependent case as well. 
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Meaning of the wave function. The Kohn-Sham system 
is noninteracting, so its total iV-particle wave function 
can be written as a single Slater determinant: 



where 



*|S(ri,...,rAr) 



1 



det{ipj} 



(23) 



The Kohn-Sham Slater determinant has only one pur- 
pose: to reproduce the exact ground-state density when 
substituted in Eq. ([6]). It is not meant to reproduce 
the exact ground-state wave function, i.e., 5'^^ 7^ in 
general. 

Meaning of the Kohn-Sham energies. The energy 
eigenvalues Sj do not have a rigorous physical meaning, 
except for the highest occupied eigenvalue. We have 

ENiN) = E{N) ~ E{N -1) = -I{N) , (24) 

i.e., the highest occupied eigenvalue of the A^-particle sys- 
tem equals minus the ionization energy of the A^-particle 
system, and 

eAr+i(A^ + 1) = E{N + 1) - E{N) = -A{N) , (25) 

i.e., the highest occupied eigenvalue of the A^-|- 1-particle 
system equals minus the electron affinity of the A^- 
particle system. 

Eigenvalue differences Ea — £», where a labels an unoc- 
cupied single-particle state and i an occupied one, should 
not be interpreted as excitation energies of the many- 
body system (although they often are). 

Asymptotic behavior of the Kohn-Sham potential. A 
neutral atom with N electrons has the nuclear poten- 
tial vo{r) = —N/r, and its Hartree potential behaves as 
vii{r) — > N/r for r — > (X). If an electron is far away in 
the outer regions of the atom, it should see the Coulomb 
potential of the remaining positive ion. This implies that 
the xc potential must behave asymptotically as 



Wxc(r) 



(26) 



for large r, for any finite system. The asymptotic behav- 
ior of z;xc(r) reflects the fact that the Kohn-Sham formal- 
ism is free of self-interaction: for a 1-electron system, the 
Hartree and xc potential cancel exactly. 

Spin-dependent formalism. In practice, the Kohn- 
Sham formalism is usually written down and applied in 
its spin-polarized form, even if the system does not have 
a net spin polarization. We then have 



V2 

— ^ + voa{r) + wh(i") + WxcCT(r) 



(27) 

where a =t, 4-. Here, the external potential voa carries 
a spin index, which could come from a static magnetic 
field, and the spin-polarized xc potential is defined as a 
functional of the spin-up and spin-down density: 



SE, 



(5rv(r) 



(28) 



(29) 



Exact exchange. The xc energy can be decomposed 
into exchange and correlation energy. The exact ex- 
change energy is given by 

TV 



Ir-r'| 



, (30) 



where the ipjaiic) are the exact Kohn-Sham orbitals. 
^oxact so-called implicit density functional. 



E. DFT in practice 

The number of applications of DFT in various ar- 
eas of science and engineering is almost impossible to 
count. Nice practical introductions from the perspectives 
of chemistry and of materials science, respectively, can 
be found in recent books and review articles [H, [ij, [l^ . 
Here, we will only make some general remarks and dis- 
cuss a couple of representative examples. 

Even though DFT is in principle exact (as we have em- 
phasized), any practical application necessarily involves 
two kinds of approximations: (i) the xc energy functional 
Exc[n], and the xc potential following from it via Eq. 
(|211) . are not exactly known and need to be approxi- 
mated; (ii) the Kohn-Sham equation (fT4|) needs to be 
solved using some computational scheme, which can in- 
troduce various types of numerical inaccuracies. 

Over the years, many approximate xc functionals 
have been proposed; some of them using physical ar- 
guments, constraints and exact conditions, others us- 
ing parametrizations combined with fitting to reference 
data. Which functional should one choose? This ques- 
tion cannot be easily answered in general but requires 
some experience. Practitioners of DFT who use popular 
software packages of quantum chemistry or solid-state 
physics often encounter daunting choices between many 
different menu options for v^c- Some functionals have 
turned out to be more popular and successful than others, 
and are typically chosen in the majority of applications. 

The xc energy of any system can be written as 



E^ 



= d 



^](r), 



(31) 



where exc['T'](r) is the xc energy density, whose depen- 
dence on the density is, in general, nonlocal: Cxc at a 
particular point r is determined by the density n(r') at 
all points in space. The goal is to approximate Cxci'^Kr). 

Much of the success of DFT can be attributed to the 
fact that a very simple approximation, the local-density 
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approximation (LDA), gives very useful results in many 
circumstances. The LDA has the following form: 



= d 



(n(r)) 



(32) 



Here, the xc energy density of a homogeneous electron 
liquid, e'^^{n) (which is a simple function of the uniform 
density n), is evaluated at the local density at point r of 
the actual inhomogeneous physical system: e'^^{n{r)) — 
e'ir{n)\- , The so defined E]l9^\ri\ is exact in the 

xt- \ ^In— n(r) l J 

limit where the system becomes uniform, and should be 
accurate when the system varies only slowly in space. 
The LDA requires e^-^ckn) as input [18]. We can write 



(33) 



where the exchange energy density can be calculated ex- 
actly using Hartree-Fock theory; the result (for the spin- 
unpolarized electron liquid) is 



1/3 



,4/3 



(34) 



This gives the following expression for the LDA exchange 
potential: 



5ri(r) 



3 /3 

4 W 



1/3 



1/3 



n(r) 



1/3 



(35) 



The correlation energy density ejf(n) is not exactly 
known, but very accurate numerical results exist from 
quantum Monte Carlo calculations. Based on these re- 
sults, parametrizations for the correlation ener gy o f the 
homogeneous electron liquid have been derived |19l - [2l| . 

The LDA generally performs very well across the 
board. It produces atomic and molecular total ground- 
state energies within 1-5% of the exact value, and yields 
molecular equilibrium distances and geometries within 
about 3%. For solids, Fermi surfaces in metals come out 
within a few percent, lattice constants of solids within 
about 2%, and vibrational frequencies and phonon ener- 
gies are obtained within a few percent as well. 

On the other hand, the LDA has several shortcom- 
ings. For instance, the LDA is not self-interaction free; 
as a consequence, the xc potential goes to zero exponen- 
tially fast and not as — 1/r [Eq. (1^5)) ]. This causes the 
Kohn-Sham energy eigenvalues to be too low in magni- 
tude in general; in particular, the highest occupied eigen- 
value Eat underestimates the ionization energy typically 
by 30-50%. The LDA does not produce any stable neg- 
ative ions, and it underestimates the band gap in solids. 
Dissociation of heteronuclear molecules in LDA produces 
ions with fractional charges. 

Overall, the LDA often gives good results in solid-state 
physics and materials science, but it is usually not accu- 
rate enough for many chemical applications. 



TABLE I. Mean absolute errors in several molecular proper- 
ties calculated for various test sets [2^ . 

Formation Ionization Equilibrium Vibrational 
enthalpy" potential*" bond length'^ frequency'' 



HF 


211.54 


1.028 


0.0249 


136.2 


LSDA 


121.85 


0.232 


0.0131 


48.9 


BLYP 


9.49 


0.286 


0.0223 


55.2 


PRE 


22.22 


0.235 


0.0159 


42.0 


B3LYP 


4.93 


0.184 


0.0104 


33.5 



" For a test set of 223 molecules (in kcal/mol). 
For a test set of 223 molecules (in eV), evaluated from the 
total-energy differences between the cation and the 
corresponding neutral, for their respective geometries. 
For a test set of 96 diatomic molecules (in A). 
For a test set of 82 diatomic molecules (in cm~^). 



The LDA can be improved by including a dependence 
not only on the local density itself, but also on gradi- 
ents of the density. This defines the so-called generalized 
gradient approximation (GGA), which has the following 
generic form: 

E^^^\n^,n^\^ J rf3^exc(nt(i-),"i(i-),Vnt(r),Vn^(r)). 

(36) 

There exist hundreds of different GGA functionals, and 
it is impossible to list all of them here. The most famous 
ones are the B88 exchange functional [23|, the LYP cor- 
relation functional [l^] (which, combined together, give 
the BLYP functional), and the PBE functional [25|. The 
exchange part of the latter has the following form: 



PBE 



1 



I +/37r2s2/3K 



(37) 



where s(r) = |Vn(r)|/2n(r)fci?(r), /ci?(r) is the local 
Fermi wavevector, and k and /3 are given parameters. 

The GGAs have been crucial in the great success story 
of DFT over the past couple of decades, due to their accu- 
racy combined with computational simplicity. However, 
improvements are still desirable. One of the most im- 
portant breakthroughs has been the development of the 
so-called hybrid functionals, which mix in a fraction of 
exact exchange: 



E. 



hybrid 



aE^ 



(I - a)E 



GGA 



E. 



GGA 



(38) 



where a is a mixing coefficient that has a value of around 
0.25. The most famous hybrid functional is B3LYP 26], 
which nowadays has become the workhorse of compu- 
tational chemistry. It should be noted that the exact 
exchange mixed in here prevents the easy construction of 
a local xc potential, so hybrid functionals are defined in 
the so-called generalized-Kohn-Sham scheme [27l[28j. 
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TABLE II. Equilibrium lattice constants of some representa- 
tive bulk solids j3Qj]. The experimental data includes a sub- 
traction of zero-point motion effects. PBEsol is a variant of 
the PRE functional [sj, and TPSS is a so-called meta-GGA 
functional, which contains a dependence on density gradients 
and on the kinetic-energy density [s^l- 

LDA PBE PBEsol TPSS Experiment 



Li 3.363 3.429 


3.428 


3.445 


3.449 


Na 4.054 4.203 


4.167 


4.240 


4.210 


Cu 3.517 3.628 


3.562 


3.575 


3.595 


Si 5.403 5.466 


5.431 


5.451 


5.416 


NaCl 5.465 5.700 


5.602 


5.703 


5.565 


MgO 4.168 4.255 


4.223 


4.237 


4.184 



Table H] gives an assessment of various approximate xc 
functionals, carried out for large sets of molecular test 
sets p2]. All xc functionals perform much better than 
Hartree-Fock. It is evident that the B3LYP functional 
gives the best overall results, with accuracies that come 
close to the requirements for predicting chemical reac- 
tions (the so-called "chemical accuracy" of 1 kcal/mol). 

In solids, hybrid functionals such as B3LYP perform 
less well, due to the fact that they do not reduce to the 
exact homogeneous electron gas limit [1^. A detailed 
assessment of the performance of modern density func- 
tionals for bulk solids was given by Czonka et al. [30j . 
Generally speaking, GGA functionals do not improve the 
lattice constants in nonmolecular solids obtained with 
LDA (which are already very good!): while LDA system- 
atically underestimates lattice constants, GGA overesti- 
mates them. Vice versa, bulk moduli and phonon fre- 
quencies are typically overestimated by LDA and under- 
estimated by GGA. This clearly affects many properties 
of solids which are volume-dependent such as their mag- 
netic behavior. Some typical results for lattice constants 
are given in Table |lTl 

A particular class of hybrid functionals, called ranqe- 
separated hybrids, has attracted much interest lately [BSj . 
The basic idea is to separate the Coulomb interaction into 
a short-range (SR) and a long-range (LR) part: 



1 



/(M|r-r'|) , l-/(M|r-r'|) 



r — r' 



r — r 



r — r 



(39) 



where the function / has the properties f{fJ-x — > 0) = 1 
and /(/ix — > oo) = 0. Common examples are f{fJ.x) — 
e~^^ and f{^ix) = erfc(/Lta;). The sep aration parameter /i 
is determined either empirically |34l - [38| or using physical 
arguments [H, [11] . The resulting range-separated hybrid 
xc functional then has the following generic form: 



P TTiSR-DFA , pLR-HF , tpDFA 



(40) 



approximation such as the LDA or GGA. The main 
strength of range-separated hybrids is that they have the 
correct (Hartree-Fock) long-range asymptotic behavior, 
and at the same time take advantage of the good short- 
range behavior of LDA or GGA. This, in turn, leads to 
a significant improvement in properties such as the po- 
larizabilities of long-chain molecules, bond dissociation, 
and, particularly importantly for TDDFT, Rydberg and 
charge-transfer excitations (see Section fVl Cp . 

This concludes our very brief survey of ground-state 
DFT. Let us now come to the dynamical case. 



III. SURVEY OF DYNAMICAL PHENOMENA 

The stationary many-body problem was defined in Sec- 
tion Solving the Schrodinger equation ([T]) allows us 
to obtain the eigenstates of an A^-particle system. The 
time-dependent Schrodinger equation is given by 

d 

« — *j(ri,...,r7v,0 = H{t)^j{vi,...,YN,t) , (41) 

where the time-dependent Hamiltonian is defined as 

H{t)^f + V{t)^W . (42) 

The time-dependent Hamiltonian has the same kinetic- 
energy and electron-electron interaction parts T and W 
as the static Hamiltonian ([2]), but it features an external 
potential operator that is explicitly time-dependent: 



V{t) 



N 



t) 



(43) 



The time-dependent Schrodinger equation (|41l) for- 
mally represents an initial-value problem. We define a 
time io our initial time (often, = 0), and we start 
with a given initial many-body wave function of the sys- 
tem, ^'(to) = '^0 (notice that this is not necessarily the 
ground state). This state is then propagated forward in 
time, describing how the system evolves under the influ- 
ence of the time-dependent potential v{v, t). In many sit- 
uations of practical interest, the time-dependent single- 
particle potential can be written as 



(44) 



i.e., the potential is static and equal to vq until time io 
when an explicitly time-dependent additional potential 
vi(t) is switched on. 

The time-dependent wave function allows us to calcu- 
late whatever observable we may be interested in. 



o{t) = {-^{t)\d\^{t)). 



(45) 



where DFA stands for any standard density-functional 



Here, 0{t) is the time-dependent expectation value of the 
Hermitian operator O corresponding to a quantum me- 
chanical observable. Two key quantities for TDDFT are 



FIG. 2. First scenario of time evolution: the external po- 
tential is static, but the system starts with a nonequilibrium 
initial state. The density then oscillates back and forth. 

the time-dependent density and current density, n(r, t) 
and j (r, t) . They can be defined via the one-particle den- 
sity operator and current density operator, 

N 

n(r)=^<5(r-r,) (46) 

1=1 

N 

j(r) = ^ [V,<5(r - r,) + <5(r - r,) V,] (47) 

1=1 

so that n{r,t) = {'i> {t)\h{r)\'9 {t)) and similar for j(r,t). 
A connection between density and current density is pro- 
vided by the continuity equation, 

i^n(r,t)--V-j(r,t). (48) 

There are many different types of quantum mechanical 
time evolution that are of practical interest. Many of 
them belong to one of the following two generic scenarios. 

First scenario. Consider a system that starts from a 
nonequilibrium initial state, and then freely evolves in 
a static potential. A simple one-dimensional example is 
illustrated in Fig. [2] at the initial time to, the density 
has an asymmetric shape which clearly does not come 
from eigenstate of the square- well potential. The density 
is then "released" and starts to oscillate back and forth, 
while the square well potential remains static {4Qi] . 

This kind of free time evolution occurs in practice when 
the system is subject to a sudden switching or a short 
"kick" at the initial time, and is then left to itself. For 
example, charge-density oscillations that are triggered in 
this way play an important role in the field of "plasmon- 
ics" 113. 

Second scenario. Consider now a system that is ini- 
tially in the ground state, and is then subject to a time- 
dependent potential that is switched on at time to. This 
is illustrated in Fig. [3] for a square well potential that is 
"shaken" by superimposing it with a time-dependent lin- 
ear potential, which again leads to an oscillating density. 

For example, this scenario takes place if an atom or 
molecule is hit by a strong laser pulse: the wave func- 
tion is driven by the external field and gets "shaken up" , 
which can then lead to ionization. 



FIG. 3. Second scenario of time evolution: the system starts 
from the ground state and evolves under a time-dependent 
external potential that is switched on at the initial time to- 



TDDFT will allow us to describe both dynamical sce- 
narios formally exactly for arbitrary many-body systems. 
To do this we will derive a dynamical version of the Kohn- 
Sham equations, which will allow us to carry out real- 
time propagations of quantum systems, starting from ar- 
bitrary initial states and under the infiuence of arbitrary 
time-dependent potentials. We will derive the formal 
framework of TDDFT in Section ITVl and we will discuss 
practical aspects and applications in Section Ivl 

Of particular importance are situations in which the 
external time-dependent potential can be considered a 
weak perturbation. Very often one is interested in the 
first-order response of the system to a perturbation, 
because many spectroscopic techniques are used this 
regime. In particular, the linear response of a material 
is directly related to its spectrum of excitations. As we 
will see in Sections IvH and Elll TDDFT in the linear- 
response regime is a very powerful approach to calculate 
excitation energies and optical spectra. In fact, this is 
where at present the majority of TDDFT applications 
are carried out. 



IV. THE FORMALISM OF TDDFT 
A. The Runge-Gross theorem 

The foundation of ground-state DFT is the Hohenberg- 
Kohn theorem, which we discussed in Section fll Bl The 
unique 1:1 correspondence between ground-state densi- 
ties and potentials makes it possible to construct den- 
sity functionals in a meaningful way, and to deter- 
mine ground-state properties in principle exactly via self- 
consistent solution of the Kohn-Sham equation. 

For the time-dependent case we would like a similar 
rigorous formal foundation. But the situation is different 
from the ground state, in two important ways. First, 
in the time-dependent case we do not have a variational 
minimum principle. Secondly, the Schrodinger equation 
(1411) is an initial- value problem, so whatever we will prove 
has to be done with a given initial state in mind. 



8 



The first to deliver an existence proof for TDDFT were 
Runge and Gross in 1984 (9|]. They proved that if two A^- 
electron systems start from the same initial state, but are 
subject to two different time-dependent potentials, their 
respective time-dependent densities will be different. 

We consider two time-dependent potentials to be dif- 
ferent if their difference is more than just a time- 
dependent constant. 



v{r,t)-v'{r,t)^c{t) 



(49) 



for t > tQ. Otherwise they would give rise to two wave 
functions that differ only by a phase factor e"*"*^*-', where 
da{t)/dt = c{t), as can easily be shown. Such purely 
time-dependent phase factors cancel out when one forms 
expectation values of operators using Eq. (^5]) . 

The Runge- Gross theorem applies to potentials that 
can be expanded in a Taylor series about the initial time: 



-v,[n](r,0 <p.(r,r) = i — ^.(r,t) ► 2]|<».(r,/)|" = «(r,0 



Density ?7(r',/') over all space and times t' <t {initial value problem) 



n(r,t') 



^0 



t 

Vs[«,%,1>o](r,0 



Time 



FIG. 4. The time-dependent Kohn-Sham equation determines 
the time-dependent density self-consistently between the ini- 
tial time to and some final time ti. The xc potential at time 
t depends on densities at times t' < t, as well as on the initial 
states of the interacting and of the Kohn-Sham system. 



v{r,t) 



OO 

E 

fe=0 



Vkjr) 
k\ 



(t-tof 



(50) 



For such potentials, the following unique 1:1 correspon- 
dence can be proven: 



unique 1:1 



n{Y,t) 



(51) 



fixed 

The proof proceeds in two steps. In the first step it is es- 
tablished that different potentials produce different cur- 
rent densities, infinitesimally later than the initial time 
^0. One then goes on to show that if the current densities 
are different, the densities must be different as well; to 
prove this, the continuity equation (|48|) is used. 

Just like in ground-state DFT, the unique 1:1 corre- 
spondence ([?T|) allows us to write the potential as a func- 
tional of the density: 



v{v,t) = v[n, *o](r,t) 



(52) 



Notice the formal dependence on the initial state. How- 
ever, this dependence goes away if the system starts from 
the ground state, i.e., Vfo = Vfgs: the Hohenberg-Kohn 
theorem then tells us that ^gsi'T-] is a functional of the 
density, and w(r, t) can be thus written as a functional of 
the density only. 

Since the potential can be written as a functional of 
the density, the time-dependent Hamiltonian becomes a 
density functional as well, and hence the time-dependent 
wave function and all observables: 

0{t) - vl/o](i)|0|*[n, v|/o](t)) = 0[n, *o](t) . (53) 



B. Time-dependent Kohn-Sham formalism 



The Kohn-Sham formalism (see Section III C|) has been 
tremendously successful in ground-state DFT. Its time- 
dependent counterpart looks very similar. The exact 



time-dependent density, n(r,t), can be calculated from 
a noninteracting system with N single-particle orbitals: 



N 



i(r,t)=^|v^j(r,i)| = 



(54) 



The orbitals (pj{r,t) satisfy the time-dependent Kohn- 
Sham equation: 



*^Vj(i-,t) 



^,(r,i), (55) 



where the time-dependent effective potential is given by 

Vs [n, *o, $o] (r, t) = w(i-, t)^v^{v, t)+v^c [n, *o, *o] (r, t) . 

(56) 

Here, v{r,t) is the time-dependent external potential, 
which we assume to have the form (j44l) . The time- 
dependent Hartree potential. 



(57) 



depends on the instantaneous time-dependent density 
only. The time-dependent xc potential formally has a 
functional dependence on the density, the initial many- 
body state ^0 of the exact interacting system, and the 
initial state of the Kohn-Sham system $o- This is 
schematically illustrated in Fig. H) 



C. Discussion: beyond Runge-Gross 

The Runge-Gross theorem in and by itself is sufficient 
to serve as the fundamental formal basis of TDDFT. 
However, there are some subtle questions that it leaves 
unanswered and some situations that are not covered by 
it. Extending the Runge-Gross theorem, or coming up 
with alternative proofs, has therefore been an area of 
significant research activity. 
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This section can be skipped by readers who may be 
less interested in the formal details of TDDFT, and more 
interested in practical aspects. 



1. v-representability and the van Leeuwen theorem 

An important question in ground-state DFT is the fol- 
lowing: given a well-behaved (i.e., continuous and not 
singular) mathematical function ri(r), with J (Prn{r) = 
N, can one always find a potential wo(r) where this 
n{r) is a ground-state density? This is known as the 
v-representability question; one distinguishes the inter- 
acting and the noninteracting v-representability problem, 
depending on whether the given density is to be repro- 
duced in the physical (interacting) or in the Kohn-Sham 
(noninteracting) system. 

Why is w-representability an important issue? If there 
exist density functions that are not v-representable (VR) , 



System 1: w(r— r') 



then the domain of the functional E„ 



would be ill de- 



fined, and one would run into formal problems in defin- 
ing functional derivatives such as in Eq. (j2ip . The v- 
representability problem in DFT is still not fully solved, 
but at least we do know that all density functions on lat- 
tice systems are VR [42] (ensemble- VR, to be precise). 
Fortunately, it turns out that the w-representability prob- 
lem in ground-state DFT can be circumvented in an ele- 
gant way with the so-called constrained search formalism 
mi [HI, which is, essentially, a clever reformulation of the 
variational minimum principle as a search over antisym- 
metric TV-particle wave functions, so that 



£;,Jn] = min*^„(«'|T + + M^l*) 



(58) 



For TDDFT the situation is different, due to a fun- 
damental difference between the ground-state problem 
and the time-dependent problem: rather than finding a 
ground state, TDDFT describes the time-propagation of 
many-body systems under the influence of external po- 
tentials. Due to the central role the external potential 
plays, the ti-representability problem [i.e., whether there 
exists a v{r,t) for every n(r,i)] seems unavoidable. 

TDDFT is not formulated on the basis of a variational 
minimum principle, since there is no quantity equivalent 
to the role of energy in time-dependent systems. In- 
stead, it is possible to formulate TDDFT via a stationary- 
action principle [45l - l47| . However, the uniqueness of 
the stationary-action point remains unproven. A rig- 
orous time-dependent version of the constrained-search 
approach does not exist, despite some attempts [48,. ,49j . 

Some progress has been made with the time-dependent 
w-representability problem for lattice systems [50|. In- 
terestingly, in this case it can happen very easily that 
perfectly well-behaved lattice densities are not VR, for 
well-understood reasons I5ll. 



The van Leeuwen theorem |52| made an important con- 
tribution towards the resolution of the w-representability 
problem in TDDFT. It makes a statement about two 




System 2: W{r-r') 



n{t) 



v'{t) uniquely 
determined 



FIG. 5. The van Leeuwen theorem states that a time- 
dependent density n{t) coming from a many-body system 
with interaction w{v — r') and potential v{v,t) can be repro- 
duced in a system with different interaction, w' (v — r') and 
potential v'{r,t). The potential v' is uniquely determined. 



many-body systems with different particle-particle inter- 
actions, w(r— r') (system 1) and w'(r— r') (system 2), see 
Fig. [5l If a time-dependent density n(r, t) is produced by 
an external potential v(r, t) in system 1 (starting from a 
given initial state), then one can uniquely construct the 
potential w'(r,t) that produces the same density in sys- 
tem 2 (the choice of initial state in system 2 is unique, 
too) . There are some restrictions on the admissible densi- 
ties: they must possess a Taylor expansion in t about the 
initial time (we denote such densities as t-TE). Below, 
we show that this assumption can be problematic. 

The van Leeuwen theorem has two important special 
cases. The first is that oi w ~ w' , i.e., the two systems 
are identical. It turns out that in this way one gains 
an alternative proof of the Runge-Gross theorem. The 
second case is w' = 0, i.e., the second system is noninter- 
acting. This establishes noninteracting u-representability 
in TDDFT, and hence provides formal justification of the 
time-dependent Kohn-Sham approach. 



2. Non-Taylor-expandable densities 

The van Leeuwen theorem shows that for f-TE densi- 
ties, one can always construct the corresponding i-TE po- 
tential for the TDKS system. However, there is a subtle 
difference between the domain of the van Leeuwen the- 
orem and the Runge-Gross theorem, as the latter only 
requires the external potentials to be t-TE [Eq. ([50]) ]. 
but not the densities. The van Leeuwen theorem does 
not apply to non-i-TE densities, which are allowed by 
the Runge-Gross theorem. Such densities are commonly 
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considered pathological and thus are not considered to 
pose any threat. However, it turns out [H, [s^ that the 
densities of most real world systems can become non-i- 
TE, including atoms, molecules, and solids! 

In the usual non-relativistic quantum mechanical de- 
scription, the nuclei and electrons interact through the 
diverging Coulomb potential, and the densities always 
have cusps at the positions of nuclei [ssj . The dynam- 
ics of the system, including the time-dependent density, 
is determined by the time evolution operator U{t,to), 
which in turn follows from the Hamiltonian H{t). In 
the presence of space-non-analytic features such as cusps, 
time-non-analyticities appear because of the kinetic en- 
ergy operator T, which is a differential operator in space. 
Thus, the time-dependent density can become non-i-TE. 

A striking example j53l] demonstrating the difference 
between the exact density and the t-TE density is shown 
in Fig. [6] (the i-TE density is defined as the result of 
applying the t-TE time-evolution operator on the initial 
state). At the initial time, a density with a cusp is pre- 
pared, and then allowed to freely evolve in time. The 
upper panel of Fig. [6] shows that the density rapidly 
becomes smooth and spreads out. By contrast, if one 
attempts to find the time evolution by using a Taylor 
expansion, the density does not move at all! 

We emphasize that although the Runge-Gross theo- 
rem is explicitly formulated for t-TE densities, the origi- 
nal proof remains valid despite the existence of non-t-TE 
densities 




[3 0.2 



sound. 



5J. Thus, the foundations of TDDFT remain 



FIG. 6. Upper panel: time-dependent density of a ID sys- 
tem with initial state ^(a;) = exp(— |a;|) propagating with no 
external potential. Lower panel: using a Taylor expansion 
in time, the initial density remains stationary. This wrong 
behavior is due to the nonanalyticity of the density. 



3. Ftxed-pomt proofs 

Recent work on the t;-representability problem and 
related questions focuses on developing so-called fixed- 
point proofs [H, [53, where the previous limitation of 
i-TE is lifted. The van Leeuwen theorem provides a way 
of constructing the time-dependent external potential for 
a given density, if the density is t-TE; if applied on non- 
i-TE densities, the constructed potential does not cor- 
respond to the exact density, but in turn reproduces the 
t-TE density [H]. The fixed-point proofs [11,113 tti^s fo- 
cus on explicitly showing the one-to-one correspondence 
between the potential and the density. The proof starts 
from the equation of motion of the density [52]: 

^^I^-W-[n{r,t)Wvir,t)]=q{r,t). (59) 

The density and the quantity q can be seen as function- 
als of the potential, and thus Eq. (j59p uniquely maps 
a potential vq to <?[wo], with the density n[i;o,^o] deter- 
mined by Wo and the initial wave function ^o- In another 
perspective, Eq. ([5^ can also be seen as a differential 
equation for the potential, when n and q are given. If this 
given density is chosen to coincide with the initial den- 
sity of the system and with its first-order time-derivative, 
and q is chosen to be q[vQ], Eq. ([5^ can be solved for 



the potential, denoted as vi. Ref. 56] proves that under 
mild restrictions, vq — vi, showing the mutual correspon- 
dence between the density and the potential. The proof 
is strengthened by recent numerical simulations [13| • The 
fixed-point proofs apply to densities confined within a fi- 
nite (but arbitrarily large) space region, and the cases 
of density cusps are included in a limiting sense. It is 
not clear as of now whether these restrictions are general 
enough for the w-representability problem. 



4- Vector potentials and time- dependent current-DFT 

TDDFT applies to electronic many-body systems in 
the presence of time-dependent scalar potentials. But 
there are important classes of time-dependent processes 
that are not included, namely, many-body systems in 
time-dependent magnetic fields or under the infiuence of 
electromagnetic waves. This is obviously a very severe 
omission, because this means that, strictly speaking, this 
precludes discussing the interaction between light and 
matter! In practice, we can often get around this restric- 
tion and treat electromagnetic fields in dipole approxima- 
tion, so that TDDFT is applicable. But in the general 
case, to deal with vector potentials of the form A{r,t) 
we need a theory that goes beyond TDDFT. 
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In general, a system can be under the influence of both 
a scalar and a vector potential, v{r,t) and A{r,t). The 
many-body Hamiltonian is then given by 



Hit) 




A(r,,t) 



+ virj,t) }+W. (60) 



The time-dependent many-body wave function associ- 
ated with H{t) determines the density n{r,t) and the 
current density j (r, t) . It is important to keep in mind 
that the current density, like any general vector field, has 
a longitudinal and a transverse component. 



(61) 



The longitudinal current density is related to the density 
via the continuity equation: 



^n(r,t) 



-V-jL(r,t), 



(62) 



but the transverse component jT(r,i) is not determined 
by n. Hence, current densities are, in general, not VR 
|58| : if j(r, t) = j£,(r,i) -|-jT(r,i) comes from a potential 
w(r,i), thenj'(r,t) = ji(r, t) -|-jy(r, t) (same longitudinal 
but different transverse component) cannot come from 
a potential w'(r,t), since this would violate the Runge- 
Gross theorem. Hence, we need the full mapping 



(u. A) o (n,j) . 



(63) 



However, this map is determined up to within a gauge 
transformation: 

z;(r,i)^w(r,t)-^A(r,t) (64) 

A(r,i) -> A(r,t) + VA(r,t) , (65) 

where A(r, t) is an arbitrary (but well-behaved) gauge 
function which vanishes at the initial time. Often, one 
chooses the gauge function in such a way that the scalar 
potential vanishes. 

Ghosh and Dhara [59l| were the first to give a formal 
proof of time-dependent current-DFT (TDCDFT). More 
recently, an alternative existence proof of TDCDFT, in 
the spirit of the van Leeuwen theorem, was provided by 
Vignale [g^l- TDCDFT on lattice systems was discussed 
by Tokatly [6lj. The time-dependent Kohn-Sham equa- 
tion in TDCDFT becomes 



this would be relativistically small. The gauge-invariant 
physical current density is given by 



j(r,t) =n(r,t)A,(r,i) + i^3[^*(r,i)V^j(r,t)] . 

(68) 



Let us summarize the key points of TDCDFT: 

1. TDCDFT overcomes formal limitations of TDDFT, 
allowing treatment of electromagnetic waves and 
general vector potentials and time-varying mag- 
netic fields. However, electromagnetic waves are 
usually treated in dipolc approximation, so one 
rarely makes use of TDCDFT in this way. 

2. The Runge-Gross theorem of TDDFT has been 
proved for finite systems, where the density van- 
ishes at infinity. However, it also works for periodic 
systems [62J , provided the external potential is also 
periodic. The Runge-Gross theorem does not ap- 
ply when a uniform homogeneous field acts on a 
periodic system. This case, however, is formally 
included in TDCDFT 0. 

3. TDCDFT can be very useful in situations that 
could, in principle, be fully described with TDDFT; 
using the current as basic variable, rather than the 
density, can make it easier to develop approxima- 
tions for dynamical xc effects [H, H^l • 



V. PRACTICAL ASPECTS 

To apply TDDFT in practice requires the following 
considerations: 

• a suitable approximation for the time-dependent xc 
potential needs to be found; 

• the time-dependent Kohn-Sham equations need to 
be solved numerically; 

• the physical observables of interest need to be ob- 
tained from the time-dependent density. 

Each of these points has its own challenges. We shall 
now address them individually, including some examples. 
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-+A,(r,t) 



+ Vsir.t) I (/3j(r,t) , 
(66) 

where the effective scalar potential, as before, is given by 
Eq. ([56]) . and the effective vector potential is 



A,(r,t) = A(r,t)+A,c(r,i) 



(67) 



Notice that the effective vector potential does not con- 
tain a Hartree-like term due to induced currents, since 



A. The time-dependent xc potential 



As we said in Section HVBl the time-dependent xc po- 
tential is formally a functional of the time-dependent den- 
sity as well as the initial states, Wxci^^, ^"01 *i'o](r, i). In 
practice, one is usually interested in situations where the 
system is initially in the ground state. If this is the case, 
things simplify considerably: thanks to the Hohenberg- 
Kohn theorem of ground-state DFT, the initial states 
become functionals of the initial (ground-state) density. 
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and the xc functional can be written as a density func- 
tional only, v^dnjiv^t). 

However, the density-dependence of the xc potential 
is complicated and nonlocal: the xc potential at space- 
time point (r, t) depends on densities at all other points in 
space and at all previous times, n{r' , t'), where t' < t (the 
potential cannot depend on densities in the future — this 
would violate the fundamental principle of causality). 

The most widely used approximation for the xc poten- 
tial is the adiabatic approximation: 



«xc(l":0 = «fcK](r)|no(r)=«(r,t), 



(69) 



where v^l, the ground-state xc potential defined in Eq. 
(j2ip . is evaluated at the instantaneous time-dependent 
density. Eq. (j69p becomes exact for an infinitely slowly 
varying system which is in its ground state for any time. 
In practice, this is of course not the case (unless one 
considers a time-dependent system which just sits there 
in its ground state, doing nothing). 

One of the most important questions in TDDFT is 
under what circumstances the adiabatic approximation 
works well. Numerical studies [65l - l67| demonstrate that 
the adiabatic approximation may break down if the sys- 
tem undergoes very rapid changes, but it turns out that 
the adiabatic approximation still works surprisingly well 
in many cases. This will be further addressed below when 
we discuss the calculation of excitation energies. 

As of today, very few applications of TDDFT have 
been carried out with nonadiabatic, explicitly memory- 
dependent xc functionals [sst - fTlj . Due to its simplic- 
ity, the overwhelming majority of time-dependent Kohn- 
Sham calculations use the adiabatic LDA (ALDA), 



,.ALDA^ 



(70) 



or any adiabatic GGA defined in a similar way, by re- 
placing the ground-state density with the instantaneous 
time-dependent density. 



B. Observables 



In Section llV Al we showed that all physical observables 
are formally functionals of the time-dependent density, 
see Eq. (|53|) . TDDFT gives, in principle, the exact time- 
dependent density n(r,t), and all quantities of interest 
must be obtained from it. Some observables are easily 
calculated in this way, but other are not. We will now 
give examples of both kinds. 



Easy observables 



quantum mechanical features such as atomic shell struc- 
ture, covalent molecular bonds, or lone pairs. Such in- 
formation can be gained from a convenient visualization 
tool known as the time-dependent electron localization 
function (TDELF) [72]. The TDELF is defined as a pos- 
itive quantity with a magnitude between zero and one: 



/ELF(r,t) 



1 



l + [D^{r,t)/D^„ir,t)Y 



The quantity 



|Vn,(r,t)|2 \jAr,tW 



8na{r,t) 2no-(r,i) 



(71) 



(72) 



is a measure of the probability of finding an electron in 
the vicinity of another electron of the same spin a at 
(r,i). Clearly, Dcr{r,t) is not an explicit density func- 
tional, but it is expressed in terms of the density, the 
current, and the orbitals via the kinetic-energy density 
rair,t) = i^fji \yipjAr,t)\\ in Eq. m is given 
by the kinetic- energy density of the homogeneous elec- 
tron liquid: 
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(67r2)3/2n5/3(j.^^) 



r.^r,t) 



(73) 



The time propagation is unitary, so the total norm is 
conserved; but to describe ionization or charge transfer 
processes, it is often of interest to obtain the number of 
electrons that escape from a given spatial region V: 



N,,,{t)=N- / dPn{v,t) 



(74) 



Here, V can be thought of as a "box" that surrounds 
the entire system (in case we wish to calculate ionization 
rates of atoms or molecules), or it could be a part of a 
larger molecule or part of a unit cell of a periodic solid. 

Another easy class of observables are moments of the 
density, such as the dipole moment: 



d(t) 



d^r rn(r, t) . 



(75) 



The dipole moment can be considered directly, i.e., in 
real time, to study the behavior of charge-density oscil- 
lations. Alternatively, it can be Fourier transformed to 
yield the dipole power spectrum |(i(a;)p or related observ- 
able quantities such as the photoabsorption cross section. 

Higher moments of the density, such as the quadrupole 
moment, can be calculated just as easily, but are less 
frequently considered. 



The easiest observable is the density itself, which shows 
how electrons move during any time-dependent process. 
This is certainly useful for visualizing molecular geome- 
tries or structural changes during chemical reactions or 
photoinduced processes, but does not reveal important 



2. Difficult observables 

Equation (|74l) gives the total number of escaped elec- 
trons, which in general can be nonintegral. For in- 
stance, if we consider an atom in a laser field, a value of 
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A^esc = 0.5 would indicate that on average half an elec- 
tron has been removed. In reality there are of course no 
"half-electrons", so we have to interpret this result in a 
probabilistic sense: it could for instance mean that there 
is 50% probability that the atom is singly ionized, and 
50% probability that it is not ionized; other scenarios, 
involving double ionization, are also possible. The prob- 
abilities to find an atom or molecule in a certain charge 
state -fm can be defined as follows I73| : 

P0(0= J d\,...J dV|*(ri,...,r^,OI' (76) 

V V 

P+i(<) = Jd\, Jd\2--.J d\N\^{r,,...,rN,t)m77) 

V V V 

and similar for all other P^™{t). Here V denotes all 
space outside the integration box V surrounding the sys- 
tem. The ion probabilities are defined in terms of the 
full many-body wave function \E'(f), which is a density 
functional according to the Runge-Gross theorem; but it 
is not possible to extract the ion probabilities P+™(i) 
directly from the density in an elementary way. 

Since the full wave function is prohibitively expensive 
to deal with, a pragmatic solution is to replace ^(t) by 
the Kohn-Sham Slater determinant $(t), in spite of the 
fact that the latter has no rigorous physical meaning. 
One then obtains the Kohn-Sham ion probabilities 

P^{t) = Niit)N2{t)...NN{t) (78) 

N 

xN^+,{t)...NN{t) (79) 
and similar for all other P+™(i), where 

N,{t) = I d'r\^,{r,t)\' . (80) 

V 

The Kohn-Sham ion probabilities are easily obtained 
from the orbitals; but apart from certain limiting cases 
[TSj , they are have no rigorous physical meaning [73, [75| . 
Here are some other examples of difficult observables: 

Photoelectron spectra. The photoelectron kinetic en- 
ergy distribution spectrum is formally defined as 

N 

P{E)dE = hm ^ \{^%\^{t))\'dE , (81) 

fc=l 

where l^*^) is a many-body eigenstate with k electron in 
the continuum and total kinetic energy E of the contin- 
uum electrons. There are approximate ways of calculat- 
ing photoelectron spectra from the density or from the 
Kohn-Sham orbitals [zilzi- 

State-to- state transition probabilities. The S-matrix 
describes the transition between two states: 

S,j = hm {^f\^{t)) , (82) 



for given initial and final many-body states '^i and 4"/. 
To get the S-matrix from the density, a cumbersome im- 
plicit read-out procedure was proposed [79j . 

Momentum distributions. Ion recoil momenta are of 
great interest in high-intense field or scattering experi- 
ments. The problem is formally similar to the problem 
of calculating ion probabilities from the density, and in 
principle requires the full wave function in momentum 
space. The Kohn-Sham momentum distributions can be 
taken as approximation, without formal justification [sof . 

Transition density matrix. The transition density ma- 
trix is a quantity that is defined in the linear response 
regime. As the name indicates, it refers to a specific 
excitation of the system (typically, a large molecular sys- 
tem), and maps the distribution and coherences of the 
excited electron and the associated hole. In particular, 
the transition density matrix is useful to visualize exci- 
tonic effects. There is no easy way to obtain it directly 
from the density; the best we can do is to construct the 
transition density matrix from Kohn-Sham orbitals [sij . 

All the above examples have in common that they are 
explicit expressions of the many-body wave function, or 
of the TV-body density matrix, but can only be implicitly 
expressed as density functionals. One can get approxi- 
mate results by replacing the full many-body wave func- 
tion with the Kohn-Sham Slater determinant ^(i), but 
there is no guarantee that this will give good results. 

C. Applications 

Real-time TDDFT has been implemented in several 
computer codes, most notably the open-source code 
octopus [Hill]. A TDDFT code must deal with two ba- 
sic numerical tasks: (i) The Kohn-Sham orbitals of the 
system, and its density, must be represented in space. 
This can be done either with a suitable basis, or on a 
spatial grid using finite-element or finite-difference dis- 
cretization schemes (octopus uses the latter), (ii) Time 
must be discretized as well, and the time-dependent 
Kohn-Sham equations are propagated forward in time, 
step by step, ensuring norm conservation. 

Let us say a few words about the time propagation. 
Suppose we know the Kohn-Sham orbitals up until some 
time T„. The orbitals at the next time step, r„+i = 
T„ -|- At, can then formally be written as 

ipj(Tn + At) = ?7(t„ + At, t„)(^j(t„) , (83) 

where t/(T„-|- At, t„) is the time evolution operator which 
propagates the orbitals one time step At forward. If At 
is sufficiently small, we can approximate U by 

{7(t„ + At, t„) « e-^^=(^"+^^/2)^- , (84) 

where Hs{Tn -\- At/2) is the time-dependent Kohn-Sham 
Hamiltonian evaluated midway between the two time 
steps (in practice, this requires a so-called predictor- 
corrector scheme [l[). The time propagation ([55)) can 
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FIG. 7. Time-dependent Kohn-Sham calculation for a CO2 
molecule. Top: time-dependent dipole moment d{t) induced 
by an initial "kick". Bottom: dipole spectrum, obtained by 
Fourier transforming d{t) (full line), compared with the spec- 
trum obtained from linear-response TDDFT (thin line). 

be numerically implemented in various ways ^84] : an ex- 
ample is the Crank-Nicholson algorithm: 

e--=- « LllSfMl , (85) 
l + iH,AT/2 

which is correct to order (At)^ and unitary (hence, the 
norm of the wave functions is conserved). This converts 
the time-dependent Kohn-Sham equations into a set of 
linear equations that can be numerically solved. 

The applications of real-time TDDFT can be roughly 
divided into two categories, related to the two scenarios 
we discussed in Section Hill 

In the first class of applications, the system is initially 
prepared in a nonequilibrium state through a sudden 
switching or a short impulsive excitation, and then al- 
lowed to propagate freely in time |85l - l87t . The initial 
perturbation is kept weak in order to avoid any non- 
linear effects, but it is spectrally broad and hence trig- 
gers a dynamical behavior of the system in which essen- 
tially the entire range of excitations participates. The 
time-dependent dipole moment d{t), Eq. ([75]). is cal- 
culated over a certain time span, and Fourier transfor- 
mation yields the optical spectrum of the system. The 
time-propagation method has certain advantages espe- 
cially for large systems [88l - [9l| and metallic clusters , 
but is less frequently used for low-lying excitations of 
smaller molecules. Below, in Section IVIl we will discuss 
an alternative way of calculating excitation energies. 

Figure [7] shows an example of such a calculation for 
the CO2 molecule. The optical absorption spectrum, 
obtained by Fourier transforming the time-dependent 
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FIG. 8. Two snapshots of the time-dependent electron local- 
ization function (TDELF) for a CO2 molecule, excited by 
a laser pulse of photon energy 20 eV and intensity 1.2 x 
10^^ W/cm^. Insets: density isosurfaces. 

dipole moment, agrees well with a spectrum that is ob- 
tained from linear-response TDDFT (we will discuss this 
approach in the following Section) . Both spectra, in turn, 
agree well with experiment [9^. 

The second class of applications is in the nonlinear 
regime, and deals with systems that are subject to strong 
excitations such as high-intensity laser pulses or collisions 
with fast, highly charged ionic projectiles. The response 
following such excitations can be highly nonlinear and 
far beyond any treatment using perturbative methods. 
Propagation of the time-dependent Kohn-Sham equa- 
tions yields the response to all orders, in principle ex- 
actly, including collective many-body effects. Quantities 
of interest include easy observables such as total ioniza- 
tion yields and high-harmonic generation spectra, and 
difficult observables such as photoelectron spectra, ion 
probabilities, or momentum distributions. 

Figure |8] shows an example. A CO2 molecule is hit with 
a very short, high-intensity laser pulse which deposits a 
large amount of excitation energy in a very short time. 
The snapshot aX t = 10.6 a.u. (1 a.u. equals 24 attosec- 
onds) shows how a packet of density flies off, and the 
remaining density is strongly distorted. The TDELF, 
Eq. ([7T|) . illustrates how the electronic orbitals have be- 
come extremely diffuse, and the bonds are essentially de- 
stroyed, which will cause the molecule to break up. 

TDDFT calculations for strong excitations have been 
carried out over the past two decades for a variety of 
atomic and molecular systems [zl [zl, HO, lo^ - foot (see 
[lOOf for a review). An intriguing question is whether 
it is possible to design the excitation (i.e., the laser in- 
tensity, pulse shape, and spectral composition) in such as 
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way that a specific control goal can be achieved. The for- 
mal framew ork of TD DFT and optimal control has been 
worked out |lOll . Il02j | , but some of the more interesting 
control goals may be difficult to achieve with standard 
(adiabatic) TDDFT approaches [T03l - [l07 |. 



VI. TDDFT AND LINEAR RESPONSE 
A. Formalism 

In many situations of practical interest, systems are 
subjected to small perturbations and hence do not de- 
viate strongly from their initial state. This happens in 
most applications of spectroscopy, where the response to 
a weak probe is used to determine the spectral proper- 
ties of a system. In this case, it is not necessary to seek a 
fully-fledged solution of the time-dependent Schrodinger 
or Kohn-Sham equations (although this would yield the 
desired information, too, as we have seen in Fig. [7|). 
Instead, one can use perturbation theory. The goal of 
linear-response theory is to directly calculate the change 
of a certain variable or observable to first order in the 
perturbation, without calculating the change of the wave 
function. For us, the most important example is the lin- 
ear density response. 

We consider the case where the system is initially in the 
ground state and a time-dependent potential is switched 
on at time tQ, see Eq. Now, however, vi{r,t) is 

treated as a small perturbation. This perturbation will 
cause some (small) time-dependent changes in the sys- 
tem, and the density will become time-dependent. We 
expand it as follows: 



n{r,t) = no(r) -I- ni(r,t) + n2(r,t) -|- . . 



(86) 



Here, uq is the ground-state density, ni is the linear den- 
sity response (the first-order change in density induced 
by the perturbation vi), n2 is the second-order density 
response (quadratic in the perturbation wi), and there 
will be higher-order terms which we have not explicitly 
indicated. If the perturbation is small, the linear density 
response dominates over all higher-order terms in the ex- 
pansion (I55|) . On the other hand, if the perturbation is 
strong, a perturbation expansion may not even converge! 
In that case it makes more sense to solve the Schrodinger 
(or Kohn-Sham) equations instead. Notice that all con- 
tributions to the density response integrate to zero, e.g., 
J (Pr ni{r,t) — 0, due to norm conservation. 

The linear density response can be formally written as 



dt' / d\'x{r,t,r\t')vi{r',t'). (87) 



ni{r,t) 



Here, x(r,r',^ — t') is the density-density response func- 
tion, defined as [l|, [l^ 

X(r, t, r', t') = -te{t - i')(*gs| [^(r, t - t'), n(r')]|vl'gs) . 



The step function 9{t — t') ensures that the response is 
causal, i.e., the response comes after the perturbation. 
Equation (|88p shows that the response function is ob- 
tained from the many-body ground state ^gg, involving 
a commutator of density operators (in interaction rep- 
resentation). Hence, via the Hohenberg-Kohn theorem, 
it is formally a functional of the ground-state density, 
x[»T-o]- Usually, one is more interested in the frequency- 
dependent response than in the real-time response: 



ni(r,w) ^ / d^r'x(r,r',cj)wi(r',w) 



(89) 



The Fourier transform of the response function (1551) can 
be written in the following form, known as the Lehmann 
representation [ll. [Tsj: 



X(r,r',w) = 



Lo — ri„ + irj 



(90) 



where the limit 77 — > 0^ is understood. Here, 

r!„ = £;„ - Eo (91) 

is the nth excitation energy of the many-body system. 
This shows explicitly that the response function has poles 
at the exact excitation energies of the system. This 
makes sense: if we apply a perturbation vi(y,u!) whose 
frequency matches one of the excitation energies, the re- 
sponse of the system is very large (we see a peak in the 
spectrum) . 

If we knew the response function x of th^ many-body 
system, calculating the density response would be easy 
and straightforward: all we have to do is evaluate ex- 
pression (j89p . From the density response, spectroscopic 
observables of interest can then be calculated. For in- 
stance, one often considers a monochromatic dipole field 
along, say, the z direction. 



vi (r, t) = £z sin(a;t) . 
The dynamic dipole polarizability follows as 

a(w) = f d^r zni{r,uj) , 



(92) 



(93) 



and the photoabsorption cross section a{uj) is given by 

(94) 



c 



In TDDFT, the linear density response can be calculated, 
in principle exactly, as the response of the noninteracting 
Kohn-Sham system to an effective perturbation [lOSj : 

ni(r,i) = j dt' j d^r'xs{r,t,r',t')vis{r',t') . (95) 
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Here, Xs(r,r',i — i') is the density-density response func- 
tion of the Kohn-Sham system. The effective perturba- 
tion is given as the sum of the real external perturbation 
plus the linearized Hartree and xc potentials: 



+ dt' dV/,e(r,t,r',t')"i(r',t')- (96) 



The so-called xc kernel is defined as the functional deriva- 
tive of the time-dependent xc potential with respect to 
the time-dependent density, evaluated at the ground- 
state density: 



/xc(r,i,r',t') 



5vy,c W(r,i) 



5n{v',t') 



(97) 



no(r) 



The effective perturbation (j96| depends on the density 
response, so the TDDFT response equation (IMI) has to 
be solved self-consistently. Again, we are usually more 
interested in the frequency-dependent response, given by 



ni(r,w) = j d3r'xs(r,r',a;)wis(r',cj) , 



(98) 



and 



(99) 



— +/xc(r,r',w) \ni{r',u}) 



The frequency-dependent xc kernel is the Fourier trans- 
form of /xc(r, t, r', t') with respect to {t — t'). 
The Kohn-Sham response function is given by 



3,k=l 



' 



(100) 

where fj and fk are occupation numbers referring to the 
configuration of the Kohn-Sham ground state (1 for oc- 
cupied and for empty Kohn-Sham orbitals), and the 
LOjk are defined as 



£k ■ 



(101) 



Thus, Xs(r,r',a;) has poles at the excitation energies 
of the noninteracting Kohn-Sham system. Naively, one 
might conclude from this that the TDDFT linear re- 
sponse must be wrong, since it contains a response func- 
tion with the wrong pole structure (we pointed out above 
that the exact response function has poles at the ex- 
act excitation energies r2„). The resolution to this ap- 
parent contradiction lies in the self-consistent nature of 
the TDDFT response equation, which "cancels out" the 
wrong poles and restores the correct poles of the many- 
body system. 

The TDDFT linear-response formalism can be gener- 
alized to a spin-dependent form. The response equation 



is then given by 

ni„{r,t) = Yl J j d^r'xsaa'ir,t,r',t')vu„ir',t') , 

(102) 

where the Kohn-Sham response function is diagonal in 
the spin index: 



X 

UJ - UJjka + IT] 



,(103) 



and ujjka = Ska — £ja- The effective perturbation is 
^;,i<,(r,w) =wi<,(r,w)+^ / c^Vi^-i— 

a' ^ 

+/xcaa'(r,r',w)|ni,,(r',L^) , (104) 
featuring the spin-dependent xc kernel /xco-o-'- 

B. How to calculate excitation energies 

The excitation energies of a many-body system are de- 
fined as the differences between the ground-state energy 
Eq and the energies of higher- lying eigenstates, £^„, see 
Eq. (|9ip. In other words, they are obtained by comparing 
the energies of stationary states. Why, then, would one 
want to use a time-dependent approach such as TDDFT? 
Isn't that unnecessarily complicated? 

It helps to think of an excitation in a different way, 
namely, as a dynamical process where the system tran- 
sitions between two eigenstates; the excitation energy 
then corresponds to a characteristic frequency, which de- 
scribes the rearrangements of probability density during 
the transition process. In other words, each excitation 
corresponds to a characteristic eigenmode of the inter- 
acting A^-electron system. 

The concept of electronic eigenmodes has a familiar 
analog in classical mechanics 10!) . A system of s coupled 
oscillators carrying out small oscillations is described by 
the homogeneous linear system of equations 

5 

^(fc„ -ri^^yOAj =0, i = l,...,s, (105) 

where the matrices kij and mij determine the potential 
and kinetic energy of the system, respectively: 



1 " 

1 * 



(106) 
(107) 
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(the Qj are generalized coordinates) . Clearly, kij and rriij 
generalize the concept of spring constant and mass of a 
simple harmonic oscillator. The solutions of Eq. (|105|) 
are obtained by finding the roots of the determinant, 



detl/cf 



= 



(108) 



The s solutions a — 1, . . . , s, are the eigenfrequen- 
cies of the system, and the associated eigenvectors Aja 
indicate the profile of the eigenmode, and can be used to 
determine the normal modes of the system. 

It turns out that calculating excitation energies with 
TDDFT is very similar to describing the small oscilla- 
tions of a classical system. Starting point is the TDDFT 
response equation, Eq. (j98l) . but without any external 
perturbation: 

n,{r,uj) = / dV'x.(r,r',w)/dV'7Hxc(r',r",w)ni(r",^) 



(109) 

where we define the combined Hartree-xc kernel as 
/Hxc(r,r',w) = |r-r'|-i + /xc(r,r',cj). Equation ((1091) 
has the trivial solution ni = for all frequencies w, but 
at certain special frequencies f2 there are also nontrivial 
solutions where the density response is finite and self- 
sustained, despite the fact that there is no external per- 
turbation. These frequencies correspond to the excita- 
tion energies of the system, and n(r, fl) is the profile of 
the associated electronic eigenmode. 

To illustrate how this works, consider the simple case 
of two electrons in a two- level system with Kohn-Sham 
orbitals (pi{r) and <y52(i"), assumed to be real. Each level 
is two-fold degenerate, and the lower level is doubly oc- 
cupied. Dropping the infinitesimal irj, the Kohn-Sham 
response function (|100p then simplifies to 



Xs{r,r',uj) 



4a;2i 

,2 _ ,.,2 



(pi(r)^2(r)^i(r')^2(r'). (HO) 



We substitute this into Eq. (|109p . and after a few simple 
manipulations we find the condition 



(111) 



where 



K{uj) = jd\j dVVi(r)(^2(r)/Hxc(r,r',w)^i(r')¥'2(r')- 

(112) 

It is a simple exercise to repeat the above example using 
the spin-dependent response formalism. Assuming that 
the ground state is not spin polarized (i.e., the spin- up 
and spin-down orbitals are the same), one finds the fol- 
lowing solutions for the eigenmodes: 



2uj2i[K„„{lo)±K„s{lo)]. 



(113) 



The plus sign represents a singlet excitation, and the mi- 
nus sign represents a triplet excitation. 

The simple examples for two-level systems are instruc- 
tive, but in practice turn out not to be quantitatively 




accurate 110l4ll2 |. The eigenmodes can be calculated, 
in pr inciple exactly, using the so-called Casida equation 
[HI: 



(114) 

where the matrix elements of A and K are given by 

Aiaa ,i' a' a' d^) = 5aa' ^aa'^aia ~^ Kiaa.i' a' a' {'-^) (115) 

X ,fiiKcaa'ir,r',uj)ip,,^,{r')ipl,^,{r') (116) 

and and a, a' run over occupied and unoccupied 
Kohn-Sham orbitals, respectively. A detailed derivation 
of Eq. piil) can be found in Ref. Q . 

If one assumes that the Kohn-Sham orbitals are real 
and that the xc kernel is frequency-independent (more 
about this assumption in section IVIDp . it is possible to 
recast the Casida equation into the following form: 



a'<T 



lacT.v a' a' 



Zi'a'a' — ■ 



(117) 



This equation can be viewed as the TDDFT counterpart 
of the eigenvalue equation (|105|) for classical small oscil- 
lations. Hence, Eq. (I117P yields the excitation energies 
and eigenmodes of the given system. 

Eq. (I114p mixes excitations and de-excitations (X and 
Y, respectively). One may simplify Eq. (|114p by setting 
the off-diagonal K matrix to zero, which decouples exci- 
tations and de-excitations. This so-called Tamm-Dancoff 
approximation (TDA) is valid if the excitation frequen- 
cies are not close to zero, which is the case for molecules, 
semiconductors, and insulators. The TDA often helps 
to compensate for deficiencies that arise because the xc 
functionals are not exactly known and have to be ap- 
proximated; the TDA can therefore be preferable over 
the full calculation (in the sense of getting qualitatively 
corre ct re sults) in certain situat ions (e.g. triplet instabil- 
ities jll4l |. conical intersections jll5j |. and excitons [ll6j ). 



C. Charge-transfer excitations 

An important class of excitations are those in which 
charge physically moves from one region (the donor) to a 
second region (the acceptor) which is spatially separated 
from the first. Such processes can occur in a wide range of 
systems, such as in complexes of two or more molecules, 
or between different functional groups within the same 
molecule. Unfortunately, the standard appro ximation s 
in TDDFT fail for charge-transfer excitations |ll7l4ll9j |. 

Consider the case where the donor and acceptor sub- 
systems are separated by a large distance R. The min- 
imum energy required to remove an electron from the 
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donor is given by the donor's ionization potential Id- 
When the electron attaches to the acceptor, some of 
that energy is regained via the acceptor's electron affinity 
Aa- Once the electron has moved from donor to accep- 
tor the two systems feel the electrostatic interaction en- 
ergy — 1/R of the induced electron-hole pair. The exact 
charge-transfer energy is therefore 



K 



(118) 



Now let us compare this with TDDFT. To make our point 
it is sufficient to consider the two-level approximation, 
Eq. ()111|) . After linearization, we obtain 

x/H.c(r,r',L.)^2(r')¥'^(r'), (HQ) 

where </?^(r) is the highest occupied donor orbital and 
f^ijc) is the lowest unoccupied acceptor orbital, which 
have exponentially vanishing overlap in the limit of large 
separation. Hence, the double integral in Eq. (|119p be- 
comes zero (assuming that the xc kernel remains finite, 
which is certainly the case for all standard approxima- 
tions), and TDDFT simply collapses to the difference 
between the bare Kohn-Sham eigenvalues. 



nr. 



(120) 



This explains why TDDFT often drastically underesti- 
mates charge-transfer excitations when co nventional xc 
functionals are used. Hybrid xc functionals |l20l . Il2l| , in 
particular the range-separated hybrids of Section lllEl of- 
fer a solution to this problem, and have been successfully 
used to des cribe cha rge-transfer excitations in a variety 
of systems [T22| - [T23 |. 



adiabatic ap proximat ion for /xc, some of the excitations 
are missing The missing excitations turn out 



D. Beyond the adiabatic approximation 

The exact excitation spectrum of a physical system is 
determined by the poles of the full response function x, 
Eq. (|90p . All of the excitation energies f2„ of the many- 
body system are, in principle, obtained by solving the 
Casida equation pi4p . But it is found that within the 
p proximat i 
|125m27| ! 

to be those that have the character of double (or mul- 
tiple) excitations, i.e., the associated many-body excited 
states, if expanded in a basis of Kohn-Sham Slater deter- 
minants, contain dominant contributions of doubly ex- 
cited configurations. 

The Kohn-Sham noninteracting response function Xs 
(jlOOl) has poles at the Kohn-Sham single excitations. 
Compared with the many-body response function 
Xs has fewer poles, since a noninteracting system can- 
not have double and multiple excitations in linear re- 
sponse. Solving the Casida equation in a finite basis and 
using the adiabatic approximation for /xc, as is done in 



practice, will not change the number of poles, but just 
shift them. To obtain double excitations, a frequency- 
dependent /xc(w) is needed which will generate addi- 
tional solutions, since the Casida equation then becomes 
a nonlinear eigenvalue problem. 

Thus, we can say the following about the adiabatic 
approximation in TDDFT: 

• The adiabatic approximation works well for those 
excitations of the physical system for which a corre- 
spondence to a single excitation in the Kohn-Sham 
system exists. The Casida equation then shifts the 
Kohn-Sham excitations towards the true single ex- 
citations. 

• The frequency dependence of /xc must kick in for 
those excitations of the physical system that are 
missing in the Kohn-Sham system, namely, double 
or multiple excitations. 

Several nonadiabatic TDDFT approaches for the de- 
scription of molecular double excitations have been ex- 
plored in the literature. One of them is known as dressed 
TDDFT [128| , where a frequency-dependent xc kernel is 
explicitly constructed within a small subspace. Other 
non adiabatic approaches are based on many-body the- 
ory [I29l4l32 |. However, none of these approaches is 
sufficiently straightforward to be part of mainstream 
TDDFT. 



E. Periodic systems and long-range behavior 



As seen from Eqs. (|115p and (I116p . the Casida equa- 
tion is expressed in th e sp ace spanned by one-particle 
Kohn-Sham transitions |l33{ . Real-space kernels are suit- 
able for calculations of finite systems such as atoms and 
molecules. For periodic systems like solids, the momen- 
tum space representation of the Hartree-xc kernel is more 
convenient. In Section IVIIDI we will use this approach 
to describe the optical properties of insulating solids. 

The real space representation of the kernel is related 
to the momentum space representation as 



/Hxc..'(r,r',..) = l E^^^"^""^- 

qeFBZ G,G' 

X /Hxc..'(q,G,G',c^)e-'('i+«')-"-', (121) 



where G, G' are reciprocal lattice vectors. With Eq. 
(I12ip . the Hartree-xc kernel in transition space, Eq. 
(|116p . becomes 



qSFBZ G,G' 

X /hxc..' (q, G, G') (a'k,, a' |e-'^(i+«')"-' |j'k,, a') 

X '5k„-ki+q,Go'5k„,-k,,+q,G^, (122) 
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with the matrix elements defined as 

(zk,a|e^('i+G)-"-|ak,a) = J Srct>l,^,{v)e'^-^+^>^ c|>a^.^,{v), 

(123) 

where k's are the Bloch wavevectors of the corresponding 
wavefunctions, and Gq, Gq can be any reciprocal lattice 
vector. The Kronecker-(5s in Eq. (|122p are a consequence 
of Bloch's theorem. 

The Hartree part of /hxc can be shown to be largely 
irrelevan t for the optical properties of insulators close to 
the gap [ l34l |: we therefore focus on the xc part in the 
following. For G = G' = (the so-called head of /xc) 
in the important limit of q — > 0, which corresponds to 
infinite range in real space, both matrix elements in Eq. 
p22p behave as 0{q^). All the local and semilocal xc 
kernels (derived from LDA and GGA in the adiabatic 
approximation) have finite values for the head. Since 
the two matrix elements in Eq. (|122p together vanish as 
0{q^), the head contribution to the sum of Eq. (|122p 
is zero for all (semi)local kernels. For these kernels, all 
changes to the Kohn-Sham spectrum come from the body 
of /xc (where G ±0^ G' ^ 0). 

Gonze et al. [l35l Il36| pointed out that the head of 
/xc has to diverge as for g — ^ to correctly de- 
scribe the polarization of periodic insulators. With the 
q^^ divergence, the head of /xc contributes in the sum 
of Eq. p22p . dominating the other parts of /xc [wings 
(G = 0,G' ^ or vice versa) and body]. Local and 
semilocal xc kernels do not have this long-range behavior, 
and there is no obvious and consistent way of modifying 
them to include the long-rangedness. 

The long-range behavior of the xc kernel is unimpor- 
tant for low-lying excitations in finite systems such as 
atoms and molecules, which means that local and semilo- 
cal xc kernels will work reasonably well. However, for 
extended and periodic systems it is crucial to have xc 
kernels with the proper long-ran ge behavior to obtain 
correct optical spectra |l34l Il37| . We will discuss this 
further in Section fVlI Dl 



VII. APPLICATIONS IN LINEAR RESPONSE 



This seemingly trivial kernel originates from many-body 
theory, where one sums up all the 'bubble' type dia- 
grams [III . Though the form is similar to time-dependent 
Hartree, TDDFT RPA is fundamentally different due to 
the use of the Kohn-Sham system. The RPA kernel has 
seen applications for molecules and is known to produce 
reasonably good results. For insulating solids, the RPA 
spectra are missing important features such as excitonic 
effects (see below). 

The proper way to obtain /xc is via Eq. (j97p : first 
approximate the time-dependent xc potential, calculate 
/xc(r, i,r',t') by taking the functional derivative, and 
then get the frequency-dependent kernel /xc(r,r',w) via 
Fourier transform. However, these steps are rarely car- 
ried out in practice, since most of the xc kernels in use are 
adiabatic kernels. Recall the adiabatic approximation for 
the xc potential, Eq. ([M]). which uses the ground-state 
functional and evaluates it at the time-dependent density. 
The adiabatic approximation for the xc kernel is 



(5no(r') (5no(r)(5no(r') ' 



(125) 



which is frequency- independent. 

An important example is the ALDA xc kernel: 



/x"c"""(r,r') = 



dv? 



5{v ~ r') , (126) 



=no(r) 



whose exchange part is explicitly given by 

/x^LDA^^^ ^ -[9^n2(r)]-i/3<5(r - r'), (127) 

and the correlation part can be obtained by using Eq. 
on any of the interpolations of ej;°^ [l3-l2l| . 
This xc kernel is not only frequency- independent, it is 
also local. One can derive adiabatic-GGA (AGGA) ker- 
nels in a similar fashion, starting from any of the standard 
GGA functionals such as those discussed in Section iHEl 
Adiabatic hybrid kernels, most notably B3LYP, are very 
widely used and have contributed much to the success of 
TDDFT in quantum chemistry. 



Linear-response TDDFT has been implemented in 
many computer codes in quantum chemistry and mate- 
rials science. In this Section we will give an overview of 
some of the most important areas of application. 



A. Standard approximations for the xc kernel 

To carry out a TDDFT calculation in the linear re- 
sponse formalism, one must know the xc kernel. The 
simplest thing to do is to use the random-phase approx- 
imation (RPA), where the xc kernel is set to zero: 

/x^/^(r,r',c.) = 0. (124) 



B. Molecular excitations 

As an example, let us consider the benzene molecule. 
Table iHll shows eight low- lying singlet and triplet excita- 
tion en ergie s of benzene, calculated with various xc func- 
tionals [l38j | . As an overall measure of the accuracy of the 
calculations, the mean absolute error (MAE) was also cal- 
culated for each functional. Based on this measure, the 
nonhybrid xc functionals (LSD, PBE, and TPSS) per- 
form at about the same level, with an MAE of 0.3-0.4 
eV. The hybrid functionals (PBEO and B3LYP) used in 
this study perform somewhat better, with an MAE rang- 
ing from 0.18 to 0.27 eV. As we will see in the following 
examples, these findings are quite typical. 
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TABLE III. Low-lying excitation energies (in eV) of the ben- 
zene molecule (CeHg) calculated with TDDFT using various 
xc functionals with the basis set 6-314— (-G(3df,3pd), and ge- 
ometry opti mize d using the respective functionals with the 
same basis ^3^. CASPT2 reference results (Ref), TDHF, 
and experimental results from Packer et al. [l39l |. The mean 
average error (MAE) for TDHF excludes the lowest (^Biu) 
triplet transition, which comes out unstable. 

LDA PBE TPSS PBEO B3LYP Ref TDHF Exp 





4.47 


3.98 


3.84 


3.68 


3.84 


3.89 




3.94 




4.82 


4.61 


4.67 


4.75 


4.72 


4.49 


4.70 


4.76 




5.33 


5.22 


5.32 


5.52 


5.41 


4.84 


5.82 


4.90 




5.05 


4.89 


4.98 


5.12 


5.07 


5.49 


5.57 


5.60 




6.07 


5.94 


6.00 


6.18 


6.05 


6.30 


5.88 


6.20 




6.12 


5.89 


5.99 


6.38 


6.11 


6.38 


6.54 


6.33 




6.70 


6.43 


6.50 


6.90 


6.62 


6.86 


6.94 


6.93 


^E2u 


6.71 


6.44 


6.50 


6.95 


6.65 


6.91 


7.11 


6.95 


MAE 


0.30 


0.37 


0.33 


0.18 


0.27 


0.09 


0.26* 
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FIG. 9. Mean absolute error for the lowest vertical excitation 
energies of a test set of 28 medium-sized organic molecules 
(103 excited states). Reproduced with permission from ACS 
from Ref. 140]. ©2009. 



Figure |9] shows the MAE for 28 xc functionals and for 
HF, obtained by calculating 103 low-lying vertical exci- 
tation ener gies for a test set of 28 medium-sized organic 
molecules |140| . compared against accurate theoretical 
benchmarks. The Kohn-Sham ground states were ob- 
tained with the same xc functionals that were used, in 
the adiabatic approximation, for the TDDFT calcula- 
tions. Identical molecular geometries were used for each 
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FIG. 10. Circular dichroism spectrum of D2-Cm, compar- 
ing TDDFT with experiment, (e: molar decadic absorption 
coefficient; R: rotatory strength; AE: excitatio n ene rgy) . Re- 
produced with permission from ACS from Ref. [144| ]. ©2002. 



xc functional. 

TDHF gives very large errors (over 1 eV), almost al- 
ways overestimating the transition energies; any TDDFT 
calculation reduces the error by at least a half. Among 
the xc functionals, we can distinguish between pure den- 
sity functionals (LDA and GGA), meta-GGAs, hybrid 
GGAs, and long-range-corrected hybrids (the first eight 
functionals in Fig. [3]). The LDA and GGAs all give an 
MAE of order 0.5 eV. Meta-GGAs (VSXC and TPSS) 
give better agreement (about 0.4 eV). But the best 
choice are clearly the hybrid GGAs (B3LYP, X3LYP, 
B98, mPWlPW91, and PBEO). In this case, the MAE is 
reduced to less than 0.25 eV. Similar findin gs w ere also 
reported in a more recent benchmark study 141 1. 

The long-range-corrected (LC) hybrids such as CAM- 
B3LYP give a slightly higher error, owing to a general 
overestimation of the transition energies. This is mainly 
due to the choice of the test set, in which charge-transfer 
excitations are not significantly represented. The advan- 
tage of long-range-corrected hybrids emerges for such ex- 
citations in larger molecules. 

As these examples illustrate, TDDFT offers an excel- 
lent compromise between computational efficiency and 
accuracy. TDDFT scales as N"^ to , depending on the 
implementation; wave- function-based methods of compa- 
rable accuracy scale at least one or two orders of magni- 
tude worse. The current limit of hig h- end wave-function- 
based methods is about 50 atoms [l40l |l42, 143]. By 
contrast, TDDFT allows the treatment of molecules con- 
taining hundreds of atoms. Examples of medium-sized 
systems are shown in Figs. [TUlandfTTl 

Figure [TU] shows the circular dichroism spectrum of a 
large chiral fullerene molecule. TDDFT was able to re- 
solve a deb ate re garding the molecular configuration of 
this system [144| . Figure [TT] shows the absorptio n spe c- 
trum of an Iridium(III) cyclometallated complex [145|. 
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FIG. 11. Calculated (blue line) and experimental (red line) 
absorption spectra of a Iridium(lII) cyclometallated complex. 
Blue vertical lines correspond to the unbroadened oscillator 
strength of the calculated singlet-singlet transitions. Repro- 
duced with permission from Elsevier from Ref. [l45l ]. ©2009. 



C. Potential-energy surfaces 

Consider a system with electrons and Nn nu- 
clei, with nuclear masses Mj and charges Zj, where 
j — 1 , . . . , Nn ■ Formally, all electrons and all nuclei 
are quantum mechanical particles, forming an interact- 
ing Ne + Nn-hody system. For instance, the H2 molecule 
depends on the coordinates of the two electrons, ri and 
r2, and on the coordinates of the two protons, Ri and 
R2: hence, it is a four-body problem. 

We denote the sets of electronic and nuclear spatial co- 
ordinates by r = {ri, . . . , tn^} and R = {Ri, . . . , RAr„}, 
respectively. The many-body eigenstates of the system 
are a function of the two sets of coordinates, 5'j(r, R), 
and obey the following many-body Schrodinger equation 



^(r,R)*»(r,R,t) = E,-^,{r,R.t) . 



(128) 



In the absence of any external potentials, the Hamilto- 
nian of the coupled electron-nuclear system is given by 



1^ 



1 



2 ^ |r, 

j.k ' J 



It 



2 ^ l-^J ^ ■^'-■1 



EE 



i=i 
Zk 



Rfc 



(129) 



As can be seen, H{t_, R) is the sum of an electronic Hamil- 
tonian containing kinetic energy and electron-electron 
interaction, Te -I- Weei a similar nuclear Hamiltonian 
T„ + Wnn, and an electron- nuclear coupling term Wen- 

The full coupled electron-nuclear many-body problem 
is too difficult to solve in general; one usually works in the 
Born-Oppenheimer (BO) approximation to obtain the 
structure of molecules and solids. The central idea of the 



BO approximation is that because of the large difference 
between the electronic and nuclear masses (the proton is 
1836 times more massive than the electron), the two sets 
of degrees of freedom are essentially decoupled. 

The BO Hamiltonian is defined as the full Hamiltonian 
(I129p minus the nuclear kinetic-energy term: 



i^Bofci)^-E^ + jE^r^ 



2 2 



2 ^ I^J ~ ■^'^l 

j^k 



E E |7~ 

j=l k=l ' ^ 



R; 



^30) 



This Hamiltonian depends parametrically on the nu- 
clear coordinates: this means that the nuclear positions 
Ri,...,RAr^ are just treated as a set of given num- 
bers, indicating a given nuclear configuration; they are 
no longer quantum mechanical operators. For each con- 
figuration one solves the Schrodinger equation 



^Bo(r,R)*j(r,R) - £;,(R)*,(r,R) . 



(131) 



The energy eigenvalues Ej{K) define the landscape of 
potential-energy surfaces, whose dimensionality depends 
on the degrees of freedom of the molecule. Thus, for a 
diatomic molecule, i?j(R) can be represented simply by a 
curve as a function of the internuclear distance, whereas 
for Nn > 3 it is a function of 3iV„ — 6 coordinates (3iV„ — 5 
for linear molecules) and should therefore more appro- 
priately be called a "hypersurface" ; the potential-energy 
surface is a 2D section through this higher-dimensional 
space. In common usage, however, the distinction be- 
tween a surface and a hypersurface is usually not made. 

The ground-state potential-energy surface i?o(R) is 
of particular interest because its minimum defines the 
molecular equilibrium position. However, excited-state 
potential-energy surfaces are important too, and play 
a crucial role in chemical reactions, photochemical pro- 
cesses, and in spectroscopy. 

All potential-energy surfaces following from Eq. (|13ip 
are called adiabatic, indicating a complete decoupling of 
electronic and nuclear degrees of freedom. The calcu- 
lation of adiabatic potential-energy surfaces is one of 
the key tasks of computational chemistry. The low- 
est potential-energy surface can be obtained exactly, 
in principle, using ground-state DFT; for excited-state 
potential-energy surfaces, forces, and vibra tional fre - 
quencies, the appropriate method is TDDFT (144 Il46| . 

Figure [12] shows the lAi manifold of the CO-stretch 
potential-energy curves of planar formaldehyde [l47| . 
These are excited states, several eV above the ground- 
state potential-energy curve (whose minimum is set at 
eV). The dashed lines are results from a multirefer- 
ence doubles CI benchmark calculation; the full lines 
were obtained with TDDFT, using the ALDA with an 
asymptotic correction. An xc functional with the correct 
asymptotics is important here because these are high- 
lying (Rydberg) excitations. 
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FIG. 12. lAi CO-stretch potential-energy curves of planar 
formaldehyde (CH2O). Full lines: TDDFT. Dashed lines: 
multireferen ce d oubles CI. Reproduced with permission from 
Wiley from [l43]. ©1998. 

A prominent feature in Fig. [T^] is the avoided crossing 
between the states labeled (tTjTt*) and {n,3py). TDDFT 
reproduces this avoided crossing qualitatively correctly, 
thanks to the configuration mixing of individual single- 
particle transitions induced by the off-diagonal matrix 
elements Kiacr,i'a'a' in the Casida equation (|114p 148]. 

The (n, 3pj,) curve is almost on top of the exact curve, 
at least for C-0 distances before the avoided crossing. 
On the other hand, the {n,3dyz) curve comes out about 
1 eV too high, primarily owing to limitations of the xc 
functional used in this calculation. 

There are many TDDFT studies in organic and inor- 
ganic photochemistry calculating excited-state potential- 
energy surfaces ,149 153]. The performance of TDDFT 
depends strongly on the xc functional used (choosing ap- 
propriate basis sets is another important factor) . Compli- 
cations can arise for potential-energy surfaces associated 
with excit ations tha t have a long-range, charge-transfer 
character jl54l . Il55| . In that case, local or semilocal xc 
functionals will fail, and one needs to use xc functionals 
with the correct long-range behavior, see Section fVl CI 

Another source of complications are situations in which 
the ground state has an intrinsically multiconfigurational 
character. This can lead to circumstances in which two 
potential-energy surfaces become degenerate and touch 
each other, which gives rise to so-called conical intersec- 
tions. The name reflects the topology in the vicinity of 
the point of degeneracy, which looks like an inverted cone 
balancing on the tip of another cone. T DDFT ha s seri - 
ous problems with conical intersections (llSl Il56l Il57l | : 



it typically produces the wrong topology in the vicinity 
of the intersection point. These difficulties have a lot 
to do with the problems of TDDFT to describe double 
excitations: an explicitly frequency-dependent xc kernel 
/xc(w) is re quire d for a proper description of conical in- 
tersections Il48ll. 



D. Optical properties of solids 

At present, the majority of applications of TDDFT are 
in the area of computational (bio)chemistry. However, 
applications in solid-state physics and materials science 
are emerging at a rapid rate. In this Section we will 
highlight some of the most important issues for TDDFT 
in solids: the band-gap problem, excitons in insulators, 
and plasmons in metals. 

1. The band gap versus the optical gap 

The fundamental band gap Eg is a key quantity that 
characterizes insulating materials. It is defined as follows: 

EgiN)^I{N)^A{N), (132) 

where I{N) and A{N) are the ionization potential and 
the electron affinity of the A''-electron system, see Eqs. 
((241) and ([25]). Hence, we obtain 

Eg{N)=eN+iiN + l)^SNiN). (133) 

It is important to note that the right-hand side of Eq. 
(|133p contains the highest occupied Kohn-Sham eigen- 
values of two different systems, namely with A'' and with 
A^ -f 1 electrons. In a macroscopic solid with 10^^ elec- 
trons, it would of course be impossible to calculate the 
band gap according to this definition. 

The band gap in the noninteracting Kohn-Sham sys- 
tem, also known as the Kohn-Sham gap, is defined as 

Eg,s{N)^eN+iiN)-SNiN). (134) 

In contrast with the interacting gap Eg , the Kohn-Sham 
gap Eg^s is simply the difference between the highest oc- 
cupied and lowest unoccupied single-particle levels in the 
same A^-particle system. This quantity is what is usually 
taken as the band gap in standard DFT band-structure 
calculations. We can relate the two gaps by 

Eg=Eg,s+A^„ (135) 

which defines Axe as a many-body correction to the 
Kohn-Sham gap. By making use of the previous rela- 
tions, we find Axe = £n+i{N + 1) — £Ar+i(A^). It turns 
out that the many-body gap correction Axe can be re- 
lated to a very fundamental property of density fu nc- 
tionals, known as derivative discontinuities 

The so-called band-gap problem of DFT reflects the 
fact that in practice Eg,s is often a poor approximation 
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Kohn-Sham gap 



Band gap 
(QP gap) 



Optical gap 



FIG. 13. Schematic illustration of the different types of 
gaps in DFT and TDDFT. The Kohn-Sham gap is defined 
as the difference of the highest occupied and lowest unoccu- 
pied Kohn-Sham eigenvalues of the A'^-electron system, see 
Eq. (|134|l . The fundamental band gap [or quasiparticle (QP) 
gap] is the Kohn-Sham gap plus the derivative discontinuity, 
see Eq. (|135p . The optical gap is the band gap minus the 
lowest exciton binding energy -E™. The Kohn-Sham gap can 
be viewed as an approximation for the optical gap. 



to -Eg, typically underestimating the exact band gap by 
as much as 50%. The reason for this is twofold: com- 
monly used approximate xc functionals (such as LDA and 
GGA) tend to underestimate the exact Kohn-Sham gap 
Eg^s, and they do not yield any discontinuity correction 
Axe • An extreme example for the second failure are Mott 
insulators, which are typically predicted to be metallic by 
DFT. This is no accident: in Mott insulators, the exact 
Kohn-Sham system is metallic (i.e., Eg^s = 0) so that 
Eg = Axe. Clearly, standard xc functionals (where Axe 
vanishes) are unfit to describe Mott insulators. 

It is important to distinguish be twee n the fundamen- 
tal band gap and the optical gap (124| . The band gap 
describes the energy that an electron must have so that, 
when added to an A^-electron system, the result is an 
N + l electron system in its ground state. The total 
charge of the system changes by —1 in this process. By 
contrast, the optical gap describes the lowest neutral ex- 
citation of an A^-electron system: here, the number of 
electrons remains unchanged. The two gaps are schemat- 
ically illustrated in Fig. [T^ltogether with the Kohn-Sham 
gap. 

The band gap of insulators can be accurately obtained 
from the so-called quasiparticle energies, which are de- 
fined as the single-particle energies of a noninteracting 
system whose one-particle Green's function is the same 
as that of the real interacting system (note how this is 
different from the definition of the Kohn-Sham system). 
In p r actice, th is is often done using the GW method 
[134 Il62l . Il63| . GW calculations are more demanding 
than DFT, but they produce band structures of solids 
that agree very well with experiment. Generalized Kohn- 
Sham schemes [13, HI] can also give good band gaps. 

While the band gap can be measured using techniques 



in which electrons are added or removed from the sys- 
tem (such as photoemission spectroscopy), the optical 
gap refers to the lowest neutral excitation. The differ- 
ence between quasiparticle band gap and optical gap is 
the lowest exciton binding energy, Eq^. Excitons can be 
viewed as bound electron-hole pairs, whose bound states 
form a Rydberg series, analogous to the hydrogen atom 
The band gap is give n by the asymptotic limit of 
the excitonic Rydberg series [l64| (at least for direct-gap 
insulators and semiconductors). 

TDDFT can be used to calculate optical spectra of ma- 
terials in principle exactly. In the case of insulators and 
semiconductors, this means that it should, in principle, 
yield the correct optical gap, the correct excitonic Ry- 
dberg series (if the material under study has one), and 
hence the correct band gap (obtained as the limit of the 
excitonic Rydberg series). We will discuss in detail in the 
following section how optical spectra of insulators and 
semiconductors are calculated with TDDFT in practice. 

As always in TDDFT, the burden rests on the xc ker- 
nel. In the case of bulk insulators, /xe needs to accom- 
plish two things: it needs to "open up" the gap (i.e., com- 
pensate the fact that the Kohn-Sham gap underestimates 
the band gap) , and it needs to produce the electron- hole 
interaction that is responsible for the forma tion of exci- 
tons. Formally, we can write this as follows jl65j |: 



f — f qp + f 



(136) 



The xc kernel is written as the sum of a quasiparticle 
part /^P (which opens up the gap) and an excitonic part 
/xc (which causes excitonic effects). The excitonic part 
turns out to be easier to approximate than the quasipar- 
ticle part (see below). In fact, no suitable approximations 
to /^P exist at present. To a large extent this is due to 
the fact that the quasiparticle part is intrinsically nona- 
diabatic [l66j : the frequency-dependence is essential to 
shift the Kohn- Sham gap, and to produce an excitonic 
Rydberg series [llj^. In view of this, one usually ignores 
the quasiparticle part of /xe and starts from a band struc- 
ture in which the gap has been corrected by other mean s 
(such as via GW, or with a simple scissor operator 167] ) . 



2. Optical spectra of semiconductors and insulators 



In the optical spectroscopy of solids, a central 
is the complex index of refraction h, defined as 



^ — ^mac 



quan tity 

in 

(137) 



where emac('^) is the macroscopic dielectric function. The 
imaginary part of emae(w) hence describes the photoab- 
sorption of a solid, as illustrated in Fig. [T3] for the 
case of silicon. To calculate the macroscopic dielectric 
function from first principles, we need to take a de- 
tour and first calculate the microscopic dielectric matrix, 
e(q, G, G',a;), where G and G' are reciprocal lattice vec- 
tors. The macroscopic dielectric function then follows as 



24 




»[eV] 



FIG. 14. Optical absorption spectrum ol bulk Si. RPA and 
TDLDA fail to reproduce the optical gap and the exc itoni c 
peak. Reproduced with permission from APS from [l37l ]. 
©2004. 



the limit [13: 



^^""^ - S e-i(q,G = 0,G' = 0,c.) 



(138) 



In turn, the inverse dielectric function of a periodic sys- 
tem can be obtained from the response function as 

e-i(q, G, G', Lo) = Sgg' + «G(q)x(q, G, G', w) , (139) 

where UG(q) = 47r/|q+Gp. In TDDFT, the fuh response 
function is expressed as 



x(q,G,G',w) = 5giG2 - Yxs{(i,Gi,G3,^^) 



G" 



G3 



X /Hxc(q, G3,G2,w) 



X.(q,G",G',c.), (140) 



GG' 



where the xc kernel in reciprocal space was defined in Eq. 
(I12ip . By calculating x on a frequency grid, one thus 
obtains the optical spectrum (including a finite broad- 
ening in order to make the spectrum smooth). The size 
of the matrices involved are determined by the number 
of k-points associated with the numerical discretization 
scheme employed. The spectral contribution from large 
G and G' elements in x typically decays rapidly, so only 
few reciprocal lattice vectors need to be considered. 

As discussed in Section IVl El the head of the xc kernel 
plays a dominant role in periodic solids. Fig. [TJ] shows 
the experimental spectrum of Si together with the calcu- 
lated spectrum of ALDA, which has a vanishing head of 
the xc matrix. Besides producing a red-shifted spectrum 
due to the band-gap problem, the ALDA spectrum lacks 
the strong excitonic peak near the gap. As expected, lo- 
cal and semilocal functionals such as the ALDA break 



down for the highly non-local excitonic effects. Big im- 
provements can be achieved by having a finite head in the 
xc kernel. We now list a few xc kernels which have the 
proper long-range behavior that is required for a finite 
head of the xc matrix. 

The long-range corrected (LRC) kernel |l36l | is a simple 
ad-hoc approximation developed mainly for studying the 
effect of the long-range behavior. It has the form 



/LRC(q,G,G',c.) = 



G|2 



(141) 



where a is a system-dependent fitting parameter. De- 
spite its simple form, LRC spectra (with proper ly chosen 
a) can be in good agreement with experiments [l37l Il69j | 
since the head contribution of the kernel usually over- 
whelms the body contributions (sometimes called local 
field effects). A simple connection of the parameter a 
with th e hig h-frequency dielectric constant has been sug- 
gested |l37l |. This xc kernel should not be confused with 
the long-range correction in ground-state DFT, where it 
means a correction term to fix the rapid dec ay o f local 
and semilocal xc potentials away from n uclei Il70l. 

The Bethe-Salpeter equation (BSE) [HI [ml is a 
many-body equation for a two-particle polarization func- 
tion (whi ch is closely related to the two-particle Green's 
function) [l73|. Today, the BSE, combined with the GW 
method, is the most accurate approach to calculating op- 
tical properties of materials. However, the scaling of the 
computational cost versus system size is not favorable; 
the use of GW-BSE has therefore been limitedto moder- 
ate system sizes, despite recent progress jl73l - ll76| . From 
the point of view of TDDFT, the BSE has been an im- 
portant guide towards the development of very accurate 
excitonic xc kernels. The idea is to construct via an 
integral equation that features the same fou r-point re - 
sponse functions that are featured in the BSE [ill [III- 
The r esul t ing xc ke rnel reproduces the results of the full 
BSE [T6i . lT78l - [l84 |. However, the computational cost is 
essentially as high as that of solving the full BSE; there- 
fore, this xc kernel has mainly served as a proof of concept 
that TDDFT is capable of producing accurate excitonic 
effects. Furthermore, the LRC xc kernel can be shown to 
emer ge fro m this BSE-based xc kernel in the long-range 
limit [18511. 

A computationally much simpler alternati ve is the re- 
cently proposed 'bootstrap' kernel |l86l . Il87| , defined as 



/,t'r'(q,G,G',c.) 



e-l(q,G,G^c. = 0)^;G,(q) 
Xo(q,G = 0,G' = 0,c<. = 0)- 



(142) 

Due to the inclusion of (q) in the numerator, the boot- 
strap kernel has the correct 0{q~'^) long-range behavior. 
The bootstrap kernel performs well for a wide range of 
solids, as illustrated in Fig. [121 and even works for the 
case of strongly bound excitons such as in solid argon or 
LiF (note that the nonintcracting response function xo 
typically contains a band-gap correction such as a scissor 
operator or GW). 
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FIG. 15. Optical absorption spectra of various bulk semi- 
conductors calculated with TDDFT using the bootstrap xc 
kerne l, Eg . (|142p . Reproduced with permission from APS 
from [ilil. ©2011. 



We also briefly mention that the VS98 meta-GGA [186 
has recently shown some promise for calculating optical 
spectra of insulators with TDDFT |i89]. 

As an alternative to obtaining optical spectra via the 
dielectric matrix, a direct calculation of excitonic binding 
energies of insulators and se miconduc tors via the Casida 
equation is also possible [53 . ll90l - IT92 |. The advantage of 
this approach is that excitonic binding energies — which 
can be in the meV range for materials such as GaAs — can 
be numerically well resolved; this is much more difficult 
to do from the dielectric function, which typically yields 
relatively low-resolution optical spectra such as in Figs. 
[T4l and [TS] It is found that the bootstrap kernel yields 
good results for strongly bound excitons , but is less ac- 
curate for the more weakly bound cases 19l| . Accurate 



triplet exciton binding energies are even more difficult to 
obtain. The development of xc kernels for excitonic ef- 
fects in solids thus remains an important task for future 
research. 

It should be noted that Eqs. ([T^ and ([T^ imply 
that the eigenvalues in the Casida equation approach 
are the poles in instead of emac, so that the absorp- 
tion peaks are not given directly. This problem is solved 
through a modification of the Hartree kernel: 



/H(q,G,G') 







G = G" = 0, 



"'^ T'^GG' otherwise. 



(143) 



|q+G| 



By using /h instead of /h in TDDFT, emac becomes [13; 
e„,ac(^) = fim [1 - i;G=o(q)x(q, G, G', w)] , (144) 

where x is the modified response function resulting from 
TDDFT with /r. Thus the Casida equation with /h 



FIG. 16. (a) Excitons arise from a coupling of single-particle 
excitations between valence and conduction band in an in- 
sulator, mediated by dynamical xc effects, (b) Particle-hole 
excitations with momentum transfer q across the Fermi sur- 
face of a simple metal. A plasmon is a coherent superposition 
of many such excitations, coupled by Coulomb interactions. 



yields eigenvalues corresponding to the peaks in the op- 
tical spectra. Since Eq. (|144p avoids the matrix inversion 
involved in Eq. (I138p . the use of /h is also a standard 
practice in the response function approach of TDDFT. 



3. Metallic systems 

The optical properties of metallic systems (bulk met- 
als or metallic nanoparticles) are strongly determined 
by the fact that they have a sea of delocalized con- 
duction electrons with a Fermi surface. Hence, their 
low-energy elementary excitation are quite different com- 
pared to systems with a gap (insulators and semiconduc- 
tors). Whereas the outstanding features of the optical 
spectra of insulators are the excitons, metallic systems 
are dominated by plasmons. 

Excitons and plasmons are observed using different ex- 
perimental techniques: excitons are seen in optical ab- 
sorption spectra (i.e., via coupling to transverse optical 
fields); on the other hand, plasmons couple to longitu- 
dinal fields, and are thus observed using electron energy 
loss spectrosc opy or in elastic light (or X-ray) scattering 
spectroscopy 

From a TDDFT perspective, both excitons and plas- 
mons are collective excitations of the many-body system. 
However, there is a big difference as to what causes the 
collective behavior in the Kohn-Sham system. Excitons 
can be viewed as a coherent superposition of a large 
number of individual particle-hole excitations between 
valence and conduct ion band, mediated via long-range 
dynamical xc effects 116j (see Fig. [TO]) . As we discussed 
in the previous subsection, it is not easy to find xc kernels 
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FIG. 17. Schematic illustration of the particle-hole continuum 
of a 3D homogeneous electron liquid, and the RPA plasmon 
dispersion. In RPA, the plasmon is undamped until it enters 
the particle-hole continuum, where it decays into incoherent 
particle- hole excitat ions (Landau damping). TDDFT gives 
very similar results [l97l ]. 
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which reproduce excitonic effects: all electron-gas based 
approximations (such as ALDA) will fail. 

On the other hand, plasmon excitations in metallic 
systems are relatively easy to capture within TDDFT. 
The reason is that plasmons can be viewed as collec- 
tive charge-density oscillations, and it is a straightfor- 
ward textbook exercise in electromagnetism to show 
that such oscillations arise already from classical elec- 
trostatic (RPA-type) interactions; many-body xc effects 
only cause relatively minor corrections (but are impor- 
tant and subtle for plasmon damping, see below). One 
thus derives the classical plasma frequency as 



Wpl 



47me^ 



TO 



(145) 



The plasmon dispersion of a homogeneous electron liquid 
can be calculated using TDDFT linear response theory, 
along similar lines as finding the zeros of the Lindhard 
dielectric function [l^ . The analytic form of the plasmon 
dispersion up to order is given by 



niq) 



1 



0, Wp/) 



(146) 

where the terms without /xc are the RPA result. For 
small q, the plasmon lies outside the particle-hole con- 
tinuum, as illustrated in Fig. [T71 As soon as the plas- 
mon dispersion enters the particle-hole continuum, it be- 
comes subject to Landau damping (decay into incoher- 
ent particle- ho le ex citations) . This damping occurs al- 
ready in RPA jl98l |. But outside the particle- hole con- 
tinuum, the only source of plasmon damping comes from 
the imaginary part of the xc kernel. The physical ori- 
gin of the low-q plasmon damping is decay into multiple 
particle-hole excitations. A frequency-independent /xc 



FIG. 18. Plasmon dispersions of bulk sodium and aluminum: 
comparison of experiment and TDDFT. Reproduced with 
permission from APS from [l9^. ©2011.) 



(such as the ALDA) has no imaginary part and hence 
leaves the plasmon undamped. 

Figure [18] shows a comparison of experimental and 
theoretical results for the plasmon dispersions of bulk 
sodium and aluminum (l96j | . The agreement is very good 
for small plasmon wavevectors, but for larger wavevec- 
tors all TDDFT approaches fai l (ev en the nonadiabatic 
xc kernel of Gross and Kohn jl08j |). Good agreement 
is achieved by a hybrid approach in which many-body 
quasiparticle lifetimes are put by hand into the response 
formahsm (TDLDA-LT). 

Plasmonic effects are found not only in bulk metals, 
but also in many types of nanostructures. TDDFT has 
been extensively used for collective excitations in metallic 
clusters and nanoparticles. In general, the results are 
excellent: plasmon peaks and line shapes for simple metal 
clusters are very well reproduced, even at the ALDA level 
(9^ Il99l |. Applications to gold and silver clusters have 
also been quite successful, and nicely demonstrate the 
evolution from atomic-like discrete spectra to plasmon 
spectra as the cluster size increases [2O0I - I202II . 

A similar picture holds for doped semiconductor nano- 
structures such as quantum wells, wires or dots. Here, 
collective excitations in the charge and spin channel have 
been well studied using TDDFT meth ods; in general, 
plasmon dispersions are well reproduced |203l | . The issue 
of plasmon damping in quantum wells has received a good 
deal of attention; in particular, intersubband plasmons in 
quantum wells have been used to te st th e Vignale-Kohn 
approximation of TDCDFT [HI [2Q|1, with consider- 
able success 
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VIII. THE FUTURE OF TDDFT 

In the final section of our overview, we attempt a fore- 
cast of the directions in which the field of TDDFT will 
be progressing. We will highlight some areas in which 
applications of TDDFT are likely to see a lot of activity 
because of their practical importance. We will also give a 
list of issues and challenges — some of them formal, some 
of them practical — which will keep the TDDFT commu- 
nity busy for years to come. 

Biological Systems. It has been said that "if the 20th 
century was the century of physi cs, the 21st century will 
be the century of biology" [209j . Without doubt, DFT 
and TDDFT methods will play a key role in the scien- 
tific effort to understand the links between structure and 
functionality in biochemistry and biology. This is due 
to the fact that DFT is the only method capable of de- 
livering ab-initio descriptions of the electronic structure 
of systems with tens of thousands of atoms; thanks to 
the development of linear-scaling methods, e ven system s 
with millions of atoms are now within reach pTo[42ll |. 

Applications of TDDFT for large biomolecu les have 
begun to emerge at a rapid rate [ol 12131 - 1219} . Many 
of these studies are concerned with the electronic and 
optical properties of DNA fragments, or the properties of 
light-harvesting complexes. Apart from the availability 
of the necessary computer power (hardware as well as 
software), there are several developments in DFT which 
facilitate this trend towards large organic systems: 

• With the range-separated hybrid functional, we 
now have the tools for describing charge-transfer 
excitations with TDDFT (see Section lyTC]) . 

• A new generation of DFT approaches for v an-der- 
Waals interactions has emerg ed m^ill, which 
allow for first-principles calculations of the struc- 
ture of sparse matter, adsorption on surfaces, and 
many other applications. 

Coupled electron-nuclear dynamics. The coupling of 
electronic and structural degrees of freedom is a decid- 
ing factor in many functionalities of biological systems. 
An example are photoinduced processes such as photoiso- 
merization. As discussed in Section fVlI C[ TDDFT gives 
access to excited-state potential energy surfaces. But 
things get really interesting when the dynamics goes be- 
yond the Born-Oppenheimer approximation, giving rise 
to effects such as structural relaxation or ultrafast laser- 
driven molecular reorganization or dissociation. In such 
situations, TDDFT can be combined with molecula r dy- 
namics, at various levels of sophistication [227l - l229| . For 
a recent review of nonadiabatic dynamics see Ref . j23Cll | . 

The most straightforward TDDFT approach for cou- 
pling electronic and nuclear dynamics is via the Ehren- 
fest approximation, which is a mixed quantum-classical 
treatment where forces on the classical ions result from a 
mean-field average over the electronic sta tes. Ehre nfest 
dynamics works well in many situations 23114234 [. but 



has its clear limitations for situations where a branch- 
ing of ionic trajectories occurs, and where the excited 
states involve multiple pathways. Such phenomena can 
be described with the so-called surface-hopping schemes 
(235 - 237j. in which multiple excited-state potential en- 
ergy surfaces can participate in the dynamics, governed 
by a stochastic hopping algorithm. 

But all of these approaches are based on classical nu- 
clear dynamics and are thus missing out on nuclear quan- 
tum effects. Important effects of nuclear dynamics such 
as interference, decoherence or tunneling are therefore 
not captured. There are already some efforts underway 
to develop approaches that combine electronic TDDFT 
with nuclear quantum dynamics 238-244]. It can be ex- 
pected that the field will continue to advance towards a 
comprehensive and practical treatment of electronic and 
nuclear degrees of freedom. This would open up a large 
area of interesting new applications of TDDFT. 

Linear and nonlinear optics in materials. In Section 
IVIIDI we discussed how linear-response TDDFT is ap- 
plied to describe optical properties of materials (insula- 
tors and metals) . It can be expected that this will remain 
a highly active area of research. Significant progress can 
be expected along several directions. 

There is a need for better xc kernels for solids. It is 
very likely that these kernels will be expressed in terms of 
occupied and unoccupied orbitals, rather than the den- 
sity. The bootstrap kernel, Eq. (|142p . is an important 
step in the right direction, but it is not so clear how it 
can be systematically improved. For instance, a spin- 
dependent generalization of the bootstrap kernel (which 
would allow a description of singlet and triplet excitons) 
is problematic (l9l| . 

A particularly hot area of research are photovoltaic 
processes in organic systems (polymers or biological light- 
harvesting complexes) [245r.248j . There is a rich va- 
riety of photophysical processes involved, such as for- 
mation and diffusion of excitons, formation of charge- 
transfer complexes, relaxation, and charge separation. At 
present, no comprehensive ab-initio picture of these pro- 
cesses exists. This represents one of the major challenges 
for TDDFT, and should soon be within reach, based on 
existing methodologies and new developments. A promis- 
ing idea is the recently proposed real-time visualization 
of exciton dynamics using the time-dependent transition 
density matrix [sif . 

In the majority of applications of TDDFT in peri- 
odic solids, the dielectric function (or related response 
properties) are calculated, which yield optical spectra or 
scattering cross sections. But there are many nonlinear 
or explicitly time-dependent processes of interest, which 
go beyond response theory and require, in principle, a 
time-dependent calculation. Real-time TDDF T calcula- 
tions for periodic solids are beginning to emerge (249l - [25^ 
to simulate hot carrier generation, dielectric breakdown, 
and coherent phonons in semiconductors and insulators. 
Such calculations, in particular if light propagation ef- 
fects are included via a coupling with Maxwell's equa- 
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tion, pose a significant computational challenge and call 
for the dev elopment of new multiscale or multidomain 
approaches |253l . l254| . 

Other developments. Let us conclude with a mixed bag 
of various formal and practical challenges and unsolved 
problems for present and future TDDFT research. 

• Nonadiabatic xc functionals. Nonadiabatic xc func- 
tionals are needed for double excitations in finite 
systems, for dissipation in extended systems, for ex- 
citon Rydberg series, for conical intersections, and 
many other important phenomena. Elect ron-g as 
based functionals [H, [1^ are of limited use |255j |: a 
connection with many-body approaches seems the 
most promising avenue towards the development of 
simple, p ractically useful nonadiabatic functionals 
129l4l32j . Another possibility co uld be via reduced 



density-matrix- functional theory j256l - l259| . 

• Open systems. TDDFT for open systems is of inter- 
est for the description of transport through nano- 
or mesoscopic systems, where a region of interest 
(e.g, a molecule) is connected to energy and parti- 
cle reservoirs via metallic leads (26d | . It is also of 
interest for treating dissipative dynamics. The cou- 
pling to a reservoir can be formally treated within 
TDDFT i n va rious ways: with a master equation 
approach j26l| , using stochastic methods [262l - l263 |. 
and by mapping the open physical^system onto a 
noninteracting closed system [265-267l|. The formal 
aspe cts a re complicated and subject of ongoing de- 
bate [268j ; practical xc functionals for open systems 
and applications beyond simple model systems can 
be expected in the future. 



• Strongly correlated systems. There has been some 
interesting recent work in which TDDFT meth- 
ods were successfully applied to the transport in 
strongly correlated model lattice systems exhibit- 



ing Coulomb blockade and the Kondo effect |269l - 
l272j | . A subtle feature of the xc potential, its deriva- 
tive discontinuity upon change of particle number 
(briefly mentioned in Section fVII Dp . turns out to 
be crucial for capturing these effects. Most of these 
studies are f or one-dim ensional Hubbard-type lat- 
tice systems |273l . |274| , but thre e -dim ensional sys- 
tems were also considered 

[M EZ^- In the fu- 
ture, work along these lines is likely to make an 
impact in the description of realistic strongly cor- 
related systems and materials, which so far have 
remained problematic for (TD)DFT. 

Extensions of the formalism. Ground-state DFT 
has long ago been extended to finite temperatures 



|277l | and to relativistic systems The corre- 

sponding TDDFT versions are not yet available, 
but would be of great interest for matter under 
extreme conditions. Finite-temperature TDDFT, 
which might include elements of nonequilibrium 
thermodynamics and time-dependent thermal en- 
sembles, could also be of interest for thermal trans- 
port and thermoelectric properties. Relativistic 
TDDFT has been used for calculating mo lecular ex - 
citation energies and response properties [27814282 1. 
and real-time Dirac -Kohn-Sham calculations have 
been explored [283| |. but formally rigorous general 
existence proofs have yet to be worked out. 
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